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5) D Claim(s) is/are allowed. 

6) |EI Claim(s) 15. 16.47-54.56-58.60 and 61 is/are rejected. 

7) E3 Claim(s) 59 is/are objected to. 

8) D Claim(s) are subject to restriction and/or election requirement. 
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DETAILED ACTION 

The Amendments filed on December 10, 2009 and September 24&25, 2009 have been 
considered. 

Petition 

The petition filed on December 10, 2009 has been granted. 

The petition filed on September 29, 2009 has been dismissed on November 23, 

2009. 

Continued Examination Under 37 CFR 1.114 

A request for continued examination under 37 CFR 1.114, including the fee set 
forth in 37 CFR 1 .1 7(e), was filed in this application after final rejection. Since this 
application is eligible for continued examination under 37 CFR 1.114, and the fee set 
forth in 37 CFR 1 .17(e) has been timely paid, the finality of the previous Office action 
has been withdrawn pursuant to 37 CFR 1.114. Applicant's submission filed on 
September 24, 2009 has been entered. 



Withdrawal of Allowability 
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The indicated allowability of claims 48, 50-54, 56-58, and 61 is withdrawn in view 
of the newly discovered reference(s) to Zhu et al. (US 2008/0294354) and Stubbs (US 
5,327,358). Rejections based on the newly cited reference(s) follow. 

Drawings 

The drawings filed on May 20, 2004 are not acceptable because: 
1/ Black shading is not acceptable, 37 CFR 1.84(m): See e.g. Figs. 10, 12, 16, 22. 
21 Character of lines, numbers, and letters are not uniformly thick and well-defined, 37 
CFR 1.84(1): See e.g. Figs. 26's, 27's, and 28-30. 

Claim Objections 

Claim 56 is objected to because of the following informalities: after "comprising" 
(line 4), "," should be -- : --; "the force" (lines 8-9) lacks antecedent basis. Appropriate 
correction is required. 

Claim Rejections - 35 USC §112 

The following is a quotation of the second paragraph of 35 U.S.C. 112: 



The specification shall conclude with one or more claims particularly pointing out and 
distinctly claiming the subject matter which the applicant regards as his invention. 
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Claims 47, 49, and 60 are rejected under 35 U.S.C. 112, second paragraph, as being 
indefinite for failing to particularly point out and distinctly claim the subject matter 
which applicant regards as the invention. 

Claims 47, 49, and 60, why is that when the "number of the stiffness parameters being 
larger than a number of system equations", "the system equations are severely 
underdetermined"? The system equations being severely underdetermined is not 
understood. 

The following is a quotation of the first paragraph of 35 U .S.C. 112: 

The specification shall contain a written description of the invention, and of the manner and process of 
making and using it, in such full, clear, concise, and exact terms as to enable any person skilled in the 
art to which it pertains, or with which it is most nearly connected, to make and use the same and shall 
set forth the best mode contemplated by the inventor of carrying out his invention. 

Claims 47, 49, and 60 are rejected under 35 U.S.C. 112, first paragraph, as 
failing to comply with the written description requirement. The claim(s) contains subject 
matter which was not described in the specification in such a way as to reasonably 
convey to one skilled in the relevant art that the inventor(s), at the time the application 
was filed, had possession of the claimed invention. The "number of stiffness 
parameters is larger than a number of system equations such that the system equations 
are severely underdetermined" is not described in the original disclosure. 



Double Patenting 
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The nonstatutory double patenting rejection is based on a judicially created 
doctrine grounded in public policy (a policy reflected in the statute) so as to prevent the 
unjustified or improper timewise extension of the "right to exclude" granted by a patent 
and to prevent possible harassment by multiple assignees. A nonstatutory 
obviousness-type double patenting rejection is appropriate where the conflicting claims 
are not identical, but at least one examined application claim is not patentably distinct 
from the reference claim(s) because the examined application claim is either anticipated 
by, or would have been obvious over, the reference claim(s). See, e.g., In re Berg, 140 
F.3d 1428, 46 USPQ2d 1226 (Fed. Cir. 1998); In re Goodman, 11 F.3d 1046, 29 
USPQ2d 2010 (Fed. Cir. 1993); In re Longi, 759 F.2d 887, 225 USPQ 645 (Fed. Cir. 
1985); In re Van Ornum, 686 F.2d 937, 214 USPQ 761 (CCPA 1982); In re Vogel, All 
F.2d 438, 164 USPQ 619 (CCPA 1970); and In re Thorington, 418 F.2d 528, 163 
USPQ 644 (CCPA 1969). 

A timely filed terminal disclaimer in compliance with 37 CFR 1 .321 (c) or 1 .321 (d) 
may be used to overcome an actual or provisional rejection based on a nonstatutory 
double patenting ground provided the conflicting application or patent either is shown to 
be commonly owned with this application, or claims an invention made as a result of 
activities undertaken within the scope of a joint research agreement. 

Effective January 1 , 1994, a registered attorney or agent of record may sign a 
terminal disclaimer. A terminal disclaimer signed by the assignee must fully comply with 
37 CFR 3.73(b). 

Claims 15, 16, 48, 50-54, 56-59, and 61 are provisionally rejected on the ground 
of nonstatutory obviousness-type double patenting as being unpatentable over claims 1 , 
3, 6, 17-20, and 36-40 of copending Application No. 12/153,348 (Zhu et al.). Although 
the conflicting claims are not identical, they are not patentably distinct from each other 
because Zhu et al. ('348) claims: 



Regarding claims 15 and 51-54, a system for determining damage information of a 
structure (claim 17), comprising: 

a sensor arranged to measure vibrations of a structure (claim 17, lines 3-4) 
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a stiffness parameter unit for receiving said vibration information, determining 
natural frequency data of said structure, and determining the stiffness parameters of 
said structure using said natural frequency data (claim 17, lines 5-8); and 

a damage information processor for receiving said stiffness parameters and 
outputting damage information comprising spatial damage information on said structure, 
said spatial damage information comprising a damage location along said lengthwise 
dimension of said structure (claims 18, 19). 

Regarding 16, said damage information processor outputs extent of damage information 
(claim 20). 

Regarding claims 48, 50, and 61 , a system for determining stiffness parameters of a 
structure (claim 1), comprising: 

a sensor arranged to measure vibrations of said structure and output vibration 
information (claim 1, lines 3-4); and 

a stiffness parameter unit for receiving said vibration information, determining 
natural frequency data of said structure, and determining the stiffness parameters of 
said structure using said natural frequency data (claim 1, lines 5-8); 

wherein said stiffness parameter unit comprises an iterative processing unit 
(claim 3) that determines said stiffness parameters using a second or higher order 
perturbation process (claim 6). 
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Regarding claim 56, a system (claim 36), comprising: 
a structure (claim 36, line 2); 

a random impact device for introducing vibrations in said structure (claim 36, 
lines 3-5), said random impact device comprising, 

a random signal generating unit for generating first and second outputs (claim 37, 
lines 3-4); 

a random impact actuator for receiving said first and second outputs (claim 37, 
lines 5-6); and 

an impact applicator coupled to said random impact actuator (claim 37, lines 7- 
8), wherein said random impact actuator drives said impact applicator such that the 
force and arrival times of said impact applicator at said structure are random (claim 37, 
lines 9-11); 

a sensor arranged to measure vibrations of said structure and output vibration 
information (claim 36, lines 6-7); and 

a stiffness parameter unit for receiving said vibration information, determining 
natural frequency data of said structure, and determining the stiffness parameters of 
said structure using said natural frequency data (claim 36, lines 8-11). 

Regarding claim 57, said random impact actuator drives said impact applicator in 
accordance with said first and second outputs (claim 38). 
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Regarding claim 58, the first and second outputs comprise independent random 
variables (claim 39). 

Regarding claim 59, the first and second outputs determine the force and arrival times, 
respectively, of the impact applicator at said structure (claim 40). 

However, regarding claim 15, Zhu et al. ('354) does not claim the structure having a 
lengthwise dimension much greater in magnitude than cross-sectional dimensions. 
However, it is obvious to determine damage of structures having a lengthwise 
dimension much greater in magnitude than cross-sectional dimensions, since these 
structures are prone to vibrating along the lengthwise direction. 

This is a provisional obviousness-type double patenting rejection because the 
conflicting claims have not in fact been patented. 

Claim Rejections - 35 USC § 102 

The following is a quotation of the appropriate paragraphs of 35 U.S.C. 102 that 
form the basis for the rejections under this section made in this Office action: 

A person shall be entitled to a patent unless - 

(b) the invention was patented or described in a printed publication in this or a foreign country or in public 
use or on sale in this country, more than one year prior to the date of application for patent in the United 
States. 
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Claims 15, 16, 48, 50-54, and 61 are rejected under 35 U.S.C. 102(b) as being 
anticipated by Stubbs (US 5,327,358). 

Regarding claims 1 5, 51 , 52, and 54, Stubbs discloses a system for determining 
damage information of a structure (Abstract, lines 1-2), comprising: 

a sensor (claim 1 , line 4) arranged to measure vibrations of a structure (claim 1 , 
lines 4-5) having a lengthwise dimension much greater in magnitude than cross- 
sectional dimensions thereof (e.g. power poles, column 1, line 16) and to output 
vibration information (measured first signal, claim 1, lines 4-5); 

a stiffness parameter unit for receiving said vibration information (column 1 , lines 
56-58; column 25, lines 31-34; 104, Fig. 5), determining natural frequency data of said 
structure (column 5, lines 8-9; column 7, lines 17-21; Table 14), and determining the 
stiffness parameters of said structure using said natural frequency data (using equation 
1, column 5, which expresses the relationship between natural frequencies and stiffness 
parameter); and 

a damage information processor (10) for receiving said stiffness parameters and 
outputting damage information comprising spatial damage information on said structure 
(column 2, line 51 - column 3, line 5), said spatial damage information comprising a 
damage location along said lengthwise dimension of said structure (column 2, lines 36- 
39). 
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Regarding claims 16, 53, and 54, Stubbs discloses said damage information processor 
outputs extent of damage information (column 2, lines 49-50). 

Regarding claims 48, 50, and 61 , Stubbs discloses a system for determining stiffness 
parameters of a structure (104, Fig. 5), comprising: 

a sensor (claim 1 , line 4) arranged to measure vibrations of said structure (claim 
1 , lines 4-5) and output vibration information (measured first signal, claim 1 , lines 4-5); 
and 

a stiffness parameter unit for receiving said vibration information (column 1 , lines 
56-58; column 25, lines 31-34; 104, Fig. 5), determining natural frequency data of said 
structure (column 5, lines 8-9; column 7, lines 17-21; Table 14), and determining the 
stiffness parameters of said structure using said natural frequency data (using equation 
1, column 5, which expresses the relationship between natural frequencies and stiffness 
parameter); 

wherein said stiffness parameter unit comprises an iterative processing unit 
(iterations, column 35, lines 16-22; column 34, lines 36-37) that determines said 
stiffness parameters using a second or higher order perturbation process (column 33, 
lines 36-54; repeating calculation of stiffness, column 34, lines 33-37). 
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Claim Rejections - 35 USC § 103 

The following is a quotation of 35 U.S.C. 1 03(a) which forms the basis for all 
obviousness rejections set forth in this Office action: 

(a) A patent may not be obtained though the invention is not identically disclosed or described as set 
forth in section 102 of this title, if the differences between the subject matter sought to be patented and 
the prior art are such that the subject matter as a whole would have been obvious at the time the 
invention was made to a person having ordinary skill in the art to which said subject matter pertains. 
Patentability shall not be negatived by the manner in which the invention was made. 

Claims 56-58 are rejected under 35 U.S.C. 103(a) as being unpatentable over Stubbs. 

Regarding claim 56, Stubbs discloses a system (Fig. 5), comprising: 

a structure (structure, Abstract, line 1 ; specimen 42); 

a random impact device (impact hammer, column 5, line 51 ) for introducing 
vibrations in said structure (column 5, lines 50-53), 

an impact applicator (impact hammer has steel tip, Google search, page 1 , 
paragraph 2) such that the force (40) and arrival times of said impact applicator at said 
structure (42) are random (column 5, lines 50-53); 

such that the force (40) and arrival times of said impact applicator at said 
structure (42) are random (column 5, lines 50-53); 

a sensor (claim 1 , line 4) arranged to measure vibrations of said structure (claim 
1 , lines 4-5) and output vibration information (measured first signal, claim 1 , lines 4-5); 
and 
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a stiffness parameter unit for receiving said vibration information (column 1 , lines 
56-58; column 25, lines 31-34; 104, Fig. 5), determining natural frequency data of said 
structure (column 5, lines 8-9; column 7, lines 17-21; Table 14), and determining the 
stiffness parameters of said structure using said natural frequency data (using equation 
1, column 5, which expresses the relationship between natural frequencies and stiffness 
parameter). 

However, Stubbs does not disclose the following claimed features: 

- Regarding claim 56, said random impact device comprising a random signal 
generating unit for generating first and second outputs; a random impact actuator for 
receiving said first and second outputs; and an impact applicator coupled to said 
random impact actuator, wherein said random impact actuator drives said impact 
applicator. 

- Regarding claim 57, said random impact actuator drives said impact applicator in 
accordance with said first and second outputs. 

- Regarding claim 58, the first and second outputs comprise independent random 
variables. 

Nevertheless, Stubbs discloses that the random impact device is a PCB board (PCB 
086B01, column 5, line 51). It would be obvious to electrically actuate the PCB impact 
device with electric signals since the device is an electrical device. 
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Therefore, it would have been obvious to a person having ordinary skill in the art at the 
time the invention was made to provide the impact device of Stubbs with electrically 
actuation means for the purpose of generating a force for inciting vibrations on a 
structure. Control of the impact device would be improved if the device is electrically 
actuated. 

Allowable Subject Matter 

Claim 59 is objected to as being dependent upon a rejected base claim, but 
would be allowable if rewritten in independent form including all of the limitations of the 
base claim and any intervening claims. 

Reasons For Allowance 

The following is an examiner's statement of reasons for allowance: 
The combination as claimed wherein a system comprising the first and second 
outputs determine the force and arrival times, respectively, of the impact applicator at 
said structure (claim 59) is not disclosed, suggested, or made obvious by the prior art 
of record. 

Any comments considered necessary by applicant must be submitted no later 
than the payment of the issue fee and, to avoid processing delays, should preferably 
accompany the issue fee. Such submissions should be clearly labeled "Comments on 
Statement of Reasons for Allowance." 
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Response to Arguments 

Applicant's arguments filed on December 10 and September 24&25, 2009 have 
been fully considered but they are not persuasive. 

With respect to the 35 USC 112, 2 nd paragraph, rejections of claims 47, 49, and 60, 
Applicants argue that "for a linear system having m equations and n unknowns, the 
system is "underdetermined" if n > m (and is "overdetermined" if m > n). Severely 
underdetermined system of linear equations include systems wherein n » m (i.e., far 
more unknowns than equations, where n represents unknowns and m represents 
equations)". 

Examiner's position is that the claims recite "the number of stiffness parameter being 
larger than a number of system equations". It is unclear whether the stiffness 
parameters are "unknown" parameters since the stiffness parameters are already 
determined by the stiffness parameter unit (see e.g. claim 47, lines 8-10). 

With respect to the 35 USC 1 12, 1 st paragraph, rejection, Applicants argue that 
"Applicant's specification discloses, inter alia, damage detection using changes of 
natural frequencies '[f]or structures such as beams and lightning masts in electric 
substations, using only the changes in the natural frequencies can relatively accurately 



Application/Control Number: 10/849,571 Page 15 

Art Unit: 2863 

detect the location(s) and extent of damage, even though the system equations are 
severely underdetermined in each iteration' (Iflf [0181]-[0182])(emphasis added) and 
discusses an example of an aluminum beam test specimen (see FIG. 12) with "severely 
underdetermined system equations (5 equations with 80 unknowns)." fl[ 
[0188])(emphasis added)". 

Examiner's position is that paragraphs [0181], [0182], and [0188] do not disclose 
comparing the system equations with the stiff parameter . Instead, paragraph [0188], 
e.g., describes comparing the system equations with the "m" unknowns (paragraph 
[0188], lines 22-24). 

Applicants further argue that "[t]he claim amendments in question were introduced in 
the Supplemental Amendment filed on December 29, 2008, and did particularly point 
out where the originally filed disclosure supported the amendments. Accordingly, the 
Examiner has failed to discharge his burden and has further failed to set forth any 
factual findings supporting the conclusory allegation of lack of written description. See, 
e.g., Purdue Pharma LP. v. Faulding Inc., 230 F.3d 1320, 1323 (Fed. Cir. 2000)(the 
written description "inquiry is a factual one and must be assessed on a case-by-case 
basis")". 

Examiner's position is that Examiner responded to the Supplemental Amendment on 
April 24, 2009 with factual findings supporting the conclusory allegation of lack of written 
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description: "Examiner's position is that paragraphs [0130] and [0188] describe 
comparing the system equations with the "m" unknowns (see paragraphs 0130, lines 
20-25; paragraph 0188, lines 22-24). However, paragraphs [0130] and [0188] do not 
disclose comparing the system equations with the stiff parameter (e.g., G_i**(0), 
paragraph 0130, line 27)" (see Office Action, filed on April 24, 2009, page 10, paragraph 
2). 

Applicant's remaining arguments have been considered but are moot in view of the new 
ground(s) of rejection. 

Contact Information 

Any inquiry concerning this communication or earlier communications from the 
examiner should be directed to Michael Nghiem whose telephone number is (571) 
272-2277. The examiner can normally be reached on M-F. 

If attempts to reach the examiner by telephone are unsuccessful, the examiner's 
supervisor, Drew Dunn can be reached on (571) 272-2312. The fax phone number for 
the organization where this application or proceeding is assigned is (571) 273-8300. 
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Information regarding the status of an application may be obtained from the 
Patent Application Information Retrieval (PAIR) system. Status information for 
published applications may be obtained from either Private PAIR or Public PAIR. 
Status information for unpublished applications is available through Private PAIR only. 
For more information about the PAIR system, see http://pair-direct.uspto.gov. Should 
you have questions on access to the Private PAIR system, contact the Electronic 
Business Center (EBC) at 866-217-9197 (toll-free). 

/Michael P. Nghiem/ 
Primary Examiner, GAU 2863 
March 1,2010 
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(54) SYSTEM AND METHOD FOR DETECTING 
STRUCTURAL DAMAGE 

(76) Inventors: Weidong Zhu, Elkridge, MD (US); 

Guangyao Xu, Baltimore, MD (US); 
Nengan Zheng, Baltimore, MD (US); 
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Related U.S. Application Data 
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20, 2003. Provisional application No. 60/512,656, 
filed on Oct. 20, 2003. 



Publication Classification 

(51) Int. CI. 7 G01H 1/14; G01M 7/08 

(52) U.S. CI 73/579; 73/12.01; 702/56 

(57) ABSTRACT 

A system and method for detecting structural damage is 
provided that utilizes a general order perturbation method- 
ology involving multiple perturbation parameters. The per- 
turbation methodology is used iteratively in conjunction 
with an optimization method to identify the stiffness param- 
eters of structures using natural frequencies and/or mode 
shape information. The stiffness parameters are then used to 
determine the location and extent of damage in a structure. 
A novel stochastic model is developed to model the random 
impact series produced manually or to generate a random 
impact series in a random impact device. The random impact 
series method or the random impact device can be used to 
excite a structure and generate vibration information used to 
obtain the stiffness parameters of the structure. The method 
or the device can also just be used for modal testing 
purposes. The random impact device is a high energy, 
random, and high signal-to-noise ratio system. 
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SYSTEM AND METHOD FOR DETECTING 
STRUCTURAL DAMAGE 

[0001] This application claims priority to U.S. Provisional 
Patent Application Ser. No. 60/471,873 filed May 20, 2003 
and U.S. Provisional Patent Application Ser. No. 60/512,656 
filed Oct. 10, 2003, which are both incorporated herein by 
reference. 

BACKGROUND OF THE INVENTION 
[0002] 1. Field of the Invention 

[0003] The invention relates to a method and apparatus for 
detecting structural damage, and, more specifically, to a 
method and apparatus for detecting structural damage using 
changes in natural frequencies and/or mode shapes. 
[0004] 2. Background of the Related Art 
[0005] Damage in a structure can be defined as a reduction 
in the structure's load bearing capability, which may result 
from a deterioration of the structure's components and 
connections. All load bearing structures continuously accu- 
mulate structural damage, and early detection, assessment 
and monitoring of this structural damage and appropriate 
removal from service is the key to avoiding catastrophic 
failures, which may otherwise result in extensive property 
damage and cost. 

[0006] A number of conventional non-destructive test 
(NDT) methods are used to inspect load bearing structures. 
Visual inspection of structural members is often unquanti- 
ii able and unreliable, especially in instances where access to 
damaged areas may be impeded or damage may be con- 
cealed by paint, rust, or other coverings. Penetrant testing 
(PT) requires that an entire surface of the structure be 
covered with a dye solution, and then inspected. PT reveals 
only surface cracks and imperfections, and can require a 
large amount of potentially hazardous dye be applied and 
disposed of. Similarly, magnetic particle testing (MT) 
requires that an entire surface of the structure be treated, can 
be applied only to ferrous materials, and detects only rela- 
tively shallow cracks. Further, due to the current required to 
generate a strong enough magnetic field to detect cracks, MT 
is not practically applied to large structures. Likewise, eddy 
current testing (ET) uses changes in the flow of eddy 
currents to detect flaws, and only works on materials that are 
electrically conductive. Ultrasonic testing (UT) uses trans- 
mission of high frequency sound waves into a material to 
detect imperfections. Results generated by all of these 
methods can be skewed due to surface conditions, and 
cannot easily isolate damage at joints and boundaries of the 
structure. Unless a general vicinity of a damage location is 
known prior to inspection, none of these methods are easily 
or practically applied to large structures which are already in 
place and operating. On the other hand, resonant inspection 
methods are not capable of determining the extent or loca- 
tion of damage, and is used only on a component rather than 
a assembled structure. None of the above NDT methods are 
easily or practically applied to large structures requiring a 
high degree of structural integrity. 

[0007] Because of these shortfalls in existing NDT meth- 
ods when inspecting relatively large structures, structural 
damage detection using changes in vibration characteristics 
has received much attention in recent years. Vibration based 
health monitoring for rotating machinery is a relatively 



mature technology, using a non-model based approach to 
provide a qualitative comparison of current data to historical 
data. However, this type of vibration based damage testing 
does not work for most structures. Rather, vibration based 
damage detection for structures is model based, comparing 
test data to analytical data from finite element models to 
detect the location(s) and extent of damage. Vibration based 
damage detection methods fall into three basic categories. 
The first of these is direct methods such as optimal matrix 
updating algorithms, which identify damage location and 
extent in a single iteration. Because of the single iteration, 
these methods are not accurate in detecting a large level of 
damage. The second category is iterative methods. The 
methodology has only been for updating modeling, which 
determines modified structural parameters iteratively by 
minimizing differences between model and test data. The 
third category includes control-based eigenstructure assign- 
ment methods, which have the similar limitation to that of 
the direct methods indicated above and are not accurate in 
detecting a large level of damage. None of these current 
vibration based methods have been incorporated into an 
iterative algorithm that can detect small to large levels of 
damage, and the vibration based approach for structures 
remains an immature technology area which is not readily 
available on a commercial basis. 

SUMMARY OF THE INVENTION 

[0008] An object of the invention is to solve at least the 
above problems and or disadvantages and to provide at least 
the advantages described hereinafter. 

[0009] Another object of the invention is to provide a 
system and method for detecting structural damage based on 
changes in natural frequencies and/or mode shapes. 
[0010] An advantage of the system and method as embod- 
ied and broadly described herein is that it can be applied to 
a large operating structure with a large number (thousands or 
more) of degrees of freedom. 

[0011] Another advantage of the system and method as 
embodied and broadly described herein is that it can accu- 
rately detect the location(s) and extent of small to large- 
levels of damage and is especially useful for detecting a 
large level of damage with severe mismatch between the 
eigenparameters of the damaged and undamaged structures. 
[0012] Another advantage of the system and method as 
embodied and broadly described herein is that it can work 
with a limited number of measured vibration modes. 

[0013] Another advantage of the system and method as 
embodied and broadly described herein is that it can use 
measurement at only a small number of locations compared 
to the degrees of freedom of the system. A modified eigen- 
vector expansion method is used to deal with the incomplete 
eigenvector measurement problem arising from experimen- 
tal measurement of a lesser number of degrees of freedom 
than that of the appropriate analytical model. 
[0014] Another advantage of the system and method as 
embodied and broadly described herein is that it can be 
applied to structures with slight nonlinearities such as open- 
ing and closing cracks. The random shaker test or the 
random impact series method can be used to average out 
slight nonlinearities and extract linearized natural frequen- 
cies and/or mode shapes of a structure. 
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[0015] Another advantage of the system and method as 
embodied and broadly described herein is that it can handle 
structures with closely spaced vibration modes, where mode 
switching can occur in the damage detection process. 

[0016] Another advantage of the system and method as 
embodied and broadly described herein is that it can handle 
different levels of measurement noise with estimation errors 
within the noise levels. 

[0017] Another advantage of the system and method as 
embodied and broadly described herein is that the damage 
detection method and the vibration testing methods such as 
the random impact series method enables damage detection 
and assessment to be automated, thus improving the reli- 
ability integrity of results. 

[0018] Another advantage of the system and method as 
embodied and broadly described herein is that damage 
detection and assessment may be automated in the field so 
that structural health can be monitored at central location 
and useful service life may be optimized. 

[0019] Another advantage of the system and method as 
embodied and broadly described herein is that the random 
impact series method enables the modal parameters such as 
natural frequencies and/or mode shapes to be measured for 
a large structure or a structure in the field when there are 
noise effects such as those arising from the wind or other 
ambient excitation. 

[0020] Additional advantages, objects, and features of the 
invention will be set forth in part in the description which 
follows and in part will become apparent to those having 
ordinary skill in the art upon examination of the following 
or may be learned from practice of the invention. The objects 
and advantages of the invention may be realized and attained 
as particularly pointed out in the appended claims. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0021] The invention will be described in detail with 
reference to the following drawings in which like reference 
numerals refer to like elements wherein: 

[0022] FIG. 1A shows a system 100 for detection of 
structural damage according to one embodiment of the 
invention; 

[0023] FIG. IB is a flowchart of an inverse algorithm for 
identifying stiffness parameters of a damaged structure from 
a select set of measured eigenparameters, in accordance with 
an embodiment of the invention; 

[0024] FIG. 2 is a flowchart of a quasi-Newton method for 
finding an optimal solution to the system equations shown in 
the flowchart of FIG. IB, in accordance with an embodi- 
ment of the invention; 

[0025] FIG. 3 is a schematic view of a serial mass-spring 

system; 

[0026] FIGS. 4A-4C illustrate estimation errors in a series 
of iterations for a low order system with a small level of 
damage; 

[0027] FIGS. 5A-5C illustrate estimation errors in a series 
of iterations for a low order system with a large level of 
damage; 



[0028] FIGS. 6A-6C illustrate estimation errors in a series 
of iterations for a large order system with a large level of 
damage; 

[0029] FIG. 7 illustrates a finite element model of a 
fixed-fixed beam; 

[0030] FIGS. 8A-8C illustrate the estimated stiffness 
parameters with complete eignenvector measurements and 
different noise levels for a ten element beam with a large 
level of damage; 

[0031] FIGS. 9A-9B illustrate stiffness parameters with 
reduced eigenvector measurements and different noise lev- 
els for a ten element beam with a large level of damage; 

[0032] FIG. 10 illustrates a modular, four bay space 
frame; 

[0033] FIG. 11 illustrates the estimated dimensionless 
stiffnesses for the damaged frame as shown in FIG. 11; 

[0034] FIG. 12 illustrates a cantilever aluminum beam test 
specimen with uniform damage in approximately 5 ele- 

[0035] FIG. 13 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of a can- 
tilever aluminum beam test specimen with the actual dam- 
age between 10 cm and 15 cm from the cantilevered end 
using a 40-element finite element model; 

[0036] FIG. 14 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of a can- 
tilevered aluminum beam test specimen with the actual 
damage between 25 cm and 30 cm from the cantilevered end 
using a 40-element finite element model; 

[0037] FIG. 15 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of an 
undamaged cantilever aluminum beam test specimen using 
a 40-element finite element model; 

[0038] FIG. 16 illustrates a cantilever aluminum beam test 

specimen with a narrow cut; 

[0039] FIG. 17 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of the test 
specimen as shown in FIG. 16 using a 45-element finite 
element model; 

[0040] FIG. 18 illustrates the estimated bending stiff- 
nesses of all the elements of the cantilever aluminum with 
simulated damage between 9 cm and 15.75 cm from the 
cantilevered end using a 40-element finite element model; 

[0041] FIG. 19 illustrates the estimated bending stiff- 
nesses of all the elements of the cantilever aluminum beam 
with simulated damage between 29.25 cm and 36 cm from 
the cantilevered end using a 40-element finite element 
model; 

[0042] FIG. 20 illustrates the estimated bending stiff- 
nesses of all the elements of the cantilever aluminum beam 
with multiple simulated damage — one at the 3 rd element and 
the other at the 20 th element; 

[0043] FIG. 21 illustrates the estimated bending stiff- 
nesses of all the elements of the cantilever aluminum beam 
with simulated uniform damage from the 16 th to the 18 th 
element; 
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[0044] FIG. 22 illustrates a 50-foot lightning mast test 
specimen with an eccentric spike at the top; 

[0045] FIG. 23 illustrates the experimental results for the 
estimated bending stiffnesses of all the elements of the 
lightning mast as shown in FIG. 22; 
[0046] FIG. 24 illustrates the estimated bending stiff- 
nesses of all the elements of the lightning mast as shown in 
FIG. 22 with simulated damage between 0.76 m and 1.125 
m from the ground using a 40-element finite element model; 
[0047] FIG. 25 illustrates the estimated bending stiff- 
nesses of all the elements of the lightning mast as shown in 
FIG. 42 with simulated damage between 6.89 m and 7.35 m 
from the ground using a 40-element finite element model; 

[0048] FIGS. 26A-26C illustrate the estimated bending 
si iffnesses of all the elements of a cantilever aluminum beam 
using the first regularization method; 

[0049] FIGS. 27A-27C illustrate the estimated bending 
stiffnesses of all the elements of a cantilever aluminum beam 
using the second regularization method; 
[0050] FIG. 28 illustrates an example of a practical appli- 
cation of the damage detection method of FIG. 1A, in 
accordance with an embodiment of the invention; 

[0051] FIG. 29 illustrates the application of a series of 
random impacts to the practical application of the damage 
detection method of FIG. 1A, in accordance with an 
embodiment of the invention; 

[0052] FIG. 30 illustrates a system for applying the series 
of random impacts shown in FIG. 29 to the practical 
application of the damage detection method of FIG. 1A, in 
accordance with an embodiment of the invention; 

[0053] FIG. 31 is a block diagram of one preferred 
embodiment of the random impact device of FIG. 30; 

[0054] FIG. 32 illustrates a lightning mast specimen; 
[0055] FIGS. 33A-33B are graphical representations of 
coherence and FRF from single and multiple impact tests of 
the mast shown in FIG. 32; 

[0056] FIG. 34 illustrates a random series of pulses with 
the same deterministic shape and random amplitudes and 
arrival times; 

[0057] FIG. 35 illustrates an average normalized shape 
function of force pulses; 

[0058] FIGS. 36-38 are graphical representations of ana- 
lytical and numerical solutions; 

[0059] FIG. 39 shows the comparison of the probability 
density function (PDF) of the experimental arrival time with 
that from uniform distribution; 

[0060] FIG. 40 shows the comparison of the experimental 
and analytical PDFs for the number of arrived pulses; and 

[0061] FIG. 41 shows the comparison of the experimental 
and analytical PDFs for the pulse amplitudes. 

DETAILED DESCRIPTION OF PREFERRED 
EMBODIMENTS 

[0062] Commonly measured modal parameters, such as 
natural frequencies and mode shapes, are functions of physi- 



cal properties of a particular structure. Therefore, changes in 
these physical properties, such as reductions in stiffness 
resulting from the onset of cracks or a loosening of a 
connection, will cause detectable changes in these modal 
parameters. Thus, if the changes in these parameters are 
indicators of damage, vibration based damage detection may 
be, simplistically, reduced to a system identification prob- 
lem. However, a number of factors have made vibration 
based damage detection difficult to implement in practice in 
the past. 

[0063] The system and method for detecting structural 
damage as embodied and broadly described herein is moti- 
vated by the observed advantages of vibration based damage 
detection over currently available technologies. It is well 
understood that this system and method may be effectively 
applied to damage detection and assessment for substantially 
all types and configurations of structures, including, but not 
limited to, simple beams, hollow tubes, trusses, frames, and 
the like. However, simply for ease of discussion, the system 
and method will first be discussed with respect to three 
examples — a mass-spring model, a beam, and a space 
frame — for conceptualization purposes. The system and 
method will later be applied lightning masts in electric 
substations. 

[0064] FIG. 1A shows a system 100 for detecting struc- 
tural damage according to one embodiment of the invention. 
System 100 includes a stiffness parameter unit 103 for 
detecting stiffness of a structure in question. Stiffness param- 
eter unit 103 may include a spectrum analyzer or any other 
device capable of receiving vibration information and pro- 
viding natural frequency and/or mode shape data. One 
example might be a spectrum analyzer capable of perform- 
ing a fast fourier transform (FFT) and with modal analysis 
software for deriving mode shape data from vibration infor- 
mation received from sensor 110 if mode shape information 
is needed. A roving sensor technique or multiple sensors 
may need to be used if mode shape information needs to be 
used. If mode shape information is desired, a sensor should 
be provided at the tip of the hammer or device used to excite 
the structure. This sensor should be configured to send data 
to the spectrum analyzer in order to obtain the frequency 
response function, which the spectrum analyzer (modal 
analysis software) uses to obtain mode shape information. 
Stiffness parameter unit 103 is configured to receive vibra- 
tion information output from the sensor 110 via sensor 
coupler 113 and input 114. Sensor coupler 113 could be a 
wire, optical fiber or even a wireless connection between 
sensor 110 and stiffness parameter unit 103. 

[0065] Stiffness parameter unit 103 may include an itera- 
tive processing unit 115 capable of determining stiffness 
parameters using a first order perturbation approach and the 
generalized inverse method. The first order perturbation 
approach can include using natural frequencies and/or mode 
shapes of the structure according to a preferred embodiment 
of the invention. Iterative processing unit 115 may be 
capable of converging to a correct result for the output of 
stiffness parameters even when the system is underdeter- 
mined. Stiffness parameter unit 103 may further include an 
outer iterative processing unit 117 and an inner (nested) 
iterative processing unit 119 which may operate using a first 
or higher order perturbation approach and a gradient or 
quasi-Newton method to be discussed below. Stiffness 
parameter unit 103 outputs stiffness parameters to damage 
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information processor 123. Damage information processor 
includes a damage extent processor 125 and a damage 
location processor 127. Damage location processor 127 
outputs the location(s) of damage and damage extent pro- 
cessor 125 outputs the extent of damage. 

[0066] System 100 using the iterative processing units 115 
or 117 and 119 is capable of determining the stiffness 
parameters and ultimately the location(s) and extent of 
damage for structures with a largest dimension less than 50 
meters, 25 meters, 10 meters, 2.5 meters and even 1.5 
meters. Examples of this are provided herein. 

[0067] The system and method may include a multiple- 
parameter, general order perturbation method, in which the 
changes in the stiffness parameters are used as the pertur- 
bation parameters. By equating the coefficients of like -order 
terms involving the same perturbation parameters in the 
normalization relations of eigenvectors and the eigenvalue 
problem, the perturbation problem solutions of all orders 
may be derived, and the sensitivities of all eigenparameters 
may be obtained. The perturbation method may be used in 
an iterative manner with an optimization method to identify 
the stiffness parameters of structures. 

[0068] Methodology 

[0069] This method presented below can simultaneously 
identify all the unknown stiffness parameters and is formu- 
lated as a damage detection problem. Since the effects of the 
changes in the inertial properties of a damaged structure are 
usually relatively small, only the changes in the stiffness 
properties due to structural damage are considered. 

[0070] Consider a N degree-of-freedom, linear, time-in- 
variant, self-adjoint system with distinct eigenvalues. The 
stiffness parameters of the undamaged structure are denoted 
by G h; (i=l, 2, . . . , m), where m is the number of the stiffness 
parameters. Structural damage is characterized by reduc- 
tions in the stiffness parameters. The estimated stiffness 
parameters of the damaged structure before each iteration 
are denoted by G ; (i=l, 2, . . . , m), and its stiffness matrix, 
which depends linearly on G ; , is denoted by K=K(G), where 
G=[G 1 , G 2 , . . . , G m ] T . Here the superscript T denotes matrix 
transpose. The eigenvalue problem of the structure with 
stiffness parameters G ; is 

K$ k =X k M<|) k (l) 

[0071] where M is the constant mass matrix, and X k =X k (G) 
and <|> k =<|> k (G) (k=l, 2, . . . , N) are the k-th th eigenvalue and 
mass-normalized eigenvector respectively. It is noted that 
X k =m k 2 , where ou k is the k-th natural frequency of the 
structure. The normalized eigenvectors of (1) satisfy the 
orthonormality relations 

(^M^&ku ($ k ) T ^ ( |)°=X k 6 ku (2) 

[0072] where liu£N and 8^ is the Kronecker delta. 
Before the first iteration, the initial stiffness parameters of 
the damaged structure are assumed to be G ; (o:) =a ; G h; (i=l, 2, 
. . . , m), where 0<a ; Sl, and the eigenvalue problem in (1) 
corresponds to that of the structure with stiffness parameters 
G ; tQ) . If there is no prior knowledge of the integrity of the 
structure, one can start the iteration from the stiffness 
parameters of the undamaged structure and set ] ; =1. Let G^ 
(i=l, 2, . . . , m) denote the stiffness parameters of the 



damaged structure. The eigenvalue problem of the damaged 
structure is 

KA d *=K*W* (3) 
[0073] where K d =K(G d ) is the stiffness matrix with G d = 
[G dl ,G d2 , . . . , G dm ] T , and X d k =X k (G d ) and <|> d k =<|> k (G d ) are 
the k-th eigenvalue and mass-normalized eigenvector 
respectively. The stiffness matrix K d is related to K through 
the Taylor expansion 




(4) 



[0074] where 8G ; =G di -Gj (i=l, 2, . . . , m) are the changes 
in the stiffness parameters, and the higher-order derivatives 
of K with respect to G ; vanish because K is assumed to be 
a linear function of G { . Based on the finite element model, 
the global stiffness matrix of a distributed structure satisfies 
(4) as its higher-order derivatives with respect to each 
element stiffness parameter vanish. 

[0075] Let the k-th eigenvalue and mass-normalized 
eigenvector of the damaged structure be related to X k and <|) k 
through 



[0076] where X (1); k , ^ (2)ij k , . . . , and X (p);j ; k are the 
coefficients of the first, second, . . . , and p-th order 
perturbations for the eigenvalue, z (1) k , z (2);j k , . . . , and z (p)ij 
. . . ; k are the coefficient vectors of the first, second, . . . , and 
p-th order perturbations for the eigenvector, and e x k and e^ 
are the residuals of order p+1. Note that the numbers in the 
parentheses in the subscripts of the coefficients and coeffi- 
cient vectors indicate the orders of the terms. By the Taylor 
expansion, one has for any p^l, 



[0077] By (7), 
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[0078] are symmetric in the p indices, i,j, . . . , and t. The 
right-hand sides of (7) are the p-th order sensitivities of the 
eigenvalues and eigenvectors with respect to the stiffness 
parameters. 



[0079] Using the normalization relations of the eigenvec- 
tors, (j) k and (|) d k , and symmetry of the coefficient vectors in 
(6), as indicated earlier, one obtains 



1 = (/ d ) T M/ d 



l + g[(/) T M4,. + (4 1)i ) r M/] 
,iC, + 

YuTu ^ 2!( * t)T * ft ™ + [ ( 4/ M 4; + ^>/ M *a>.] +2!(^„/M*'}aG 1 aG J - + 



2![(*fiu/Afcfi)i + (4>./ M 4./ + + 3 !(z, 3)i:/1 ) W ]SG,SGj6G, + 



2!(p - 2)![(4 2W ) T M4_ 211 ... J( + (4r w l'*'.4,. 2)j » + " + (4»/ M 4-2)!,- r] + - + 
(P - 2)!2![(^ 2y ... a ) T Aftf 2)ff + (4_ 21j ... J r Mzf 21i , + - + (4-2),,.../^,] + 
(p - 1)![(^ 1W ...„) T M4 V + (4- lvl ... jM^j + - + (4_ 1)i ,-.../A/4 1 ,] + 

P'[(4« ...afMiffiGiSGj ■■■ SG S SG, + ■■■ 
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[0080] where the superscript T denotes transpose. For any 
p-th (p = l) order term in the last expression of (8), the 



4-., 



[0081] is defined by 

R} j ... t =X l \X 2 \---X a \, 



-continued 

(42)22) 7 M42)12 + (42)22f ^2)12 + 

Hi&WlfM^ + (4, 12 /a/4 2 + 

(4l22) 7 'M4.|2+(4) 1 22) 7 'M4 1 .2] + 

4U;i hl222 : r MJ, k }0G l SGi = 
{4(0* ) T he tw222 + [(4,,, ) r M4c2: + 3(4, )2 > T Mz* 3)122 ] + 

4(42)12) r M42)22 + [(4)222) T ^4l)l + 

3! 4,122 l?M 4l 12 I + 4 <4l 11222 > 7 M 0* )<SG| iSG. 



[0082] where X 1; X 2 , . . . , and ^ (liaSp) are the 
numbers of the first, second, . . . , and last distinct indices 
within indices, i, j, . . . , and t, with X 1 +X 2 + . . . +X a =p. For 
instance, R 1234 -l!l!l!l!-l with a=4, and R 112223 =2!3!l!= 
2!3! with a=3. Some explanations of the general p-th order 
term in the last expansion in (8) are in order. It includes p+1 
types of terms ordered from the beginning to the end of the 
expression within the braces, with each type of terms except 
the first and last onus enclosed in square brackets. The c-th 
(l^c^p+1) type of terms is obtained by multiplying a 
(c-l)-th order term in the expansion of (<|> d k ) T by a (p-c+1)- 
th order term in the expansion of M(|> d k . For the c-th type of 
term, where 2^cSp, the p indices, i, j, . . . , and t are 
distributed in the subscripts of the two vectors pre- and 
post-multiplying M, whose numbers of indices in the sub- 
scripts are c-1 and p-c+1 respectively. For the first and last 
((p+l)-th) types of terms, all the p indices lie in the sub- 
scripts of the vectors post- and pre-multiplying M, respec- 
tively. The number of terms within each set of square 
brackets equals the number of different combinations of 
indices in the subscripts of the vectors pre- and post- 
multiplying M. When all the p indices, i,j, . . . , and t, have 
distinct values, due to symmetry of these vectors in their 
indices, terms of the c-th (l^c^p+1) type, involving dif- 
ferent permutations of the same indices in the subscripts of 
the vectors, are equal and combined, resulting in the scaling 
factor (c-l)!(p-c+l)! in front of the square brackets. Con- 
sequently, the indices in the second through p-th summations 
range from the values of their preceding summation indices 
to m. When any of the p indices, i, j, . . . , and t, have equal 
values, the corresponding terms in each type are given by 
those in the previous case divided by 



[0083] For instance, the 4 th order term of the form 
5G 1 6G 2 3 in the expansion of (i|) d k ) T M(|> d k can be obtained 
from the expression for the p-th order term in (8): 



[0084] where p=4 and the four indices involved, i, j, 1, and 
o, are i=l and j=l=o=2. 

[0085] liquating the coefficients of the 60, (i=l, 2, . . . , m) 
terms in (8) and using symmetry of the mass matrix yields 
(<ft T M Z(1 /=0 (10) 
[0086] Equating the coefficients of the SGjSGj terms and 
using symmetry of M and z (2)ij k yields 



[0087] for all i, j=l,2, . 
procedure, one obtains 



, m . Following a similar 



-^[tfivf M $)jl + (4l)/ M 42)i; + (4i)/m42 W 



[0088] for i,j,l=l, 2, . . . , m. Equating the coefficients of 
the bG fiGj . . . 50,50, terms of p-th order in (8) yield 



(/) r M4 )ff ...,= 1 

"^{(P - Dllftfi/Mi^iur...- + tfij/jftUw-- + 
2!(p - 2)![(zf 2) ,j.) r A/; ( ,,_ 21 , ... s , + (^ 2v f M^ 2)j ... „ + 

■■■+(W' M 4-2)i/...,]+-" + 

(p-2)-2-[(4_ 2)1 ... K )''M4 2e + (4_ 2)j ,.. j ,)'m4 2) , + 
(p-l)![(4^i W ... Jm^ + ^W- J'M^j + 



T ^{4!(/fM4 1222 + 3-[(4, ) ,) T M43 ) 222 + (4.2) 7 M43 ) i : 

(4ll2) r M4 3 )122 + (4l)2) 7M 43)122] + 

2!2![(4 )12 ) 7 M4 2)22 + + 



[0089] for i,j, . . . ,t=l,2, . . . , m. The p-1 types of terms, 
enclosed in the p-1 sets of square brackets on the right-hand 
side of (13), are ordered from the beginning to the end of the 
expression within the braces, and their structures are readily 
observed. When p is odd, by symmetry of 
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Mand 4„ y ...„the y - r/.(l < y < ?-J-} 



[0090] type of terms is identical to the (p-y)-th type, and 
the two types of terms can be combined. Similarly, when p 



[0097] is the coefficient of the u-th term and the number ir 
the parentheses in its subscript corresponds to that of 



[0091] type of terms equals and can be combined with the 
(p-y)-th type of terms. In this case, however, there is a 
separate type, the 



[0098] Pre-multiplying (15) by (<|) k ) r and using (1), (2), 

and (16) yields 



[0092] of terms in the middle of the expression. 
[0093] Substituting (4)-(6) into (3) yields 



[0099] Substituting (16) into (10) and using (2) yields 



ggg^<w*---}= 

|a» + g A< V SG: + g g ^jSCSGj + 
i=i j=i (=i J 

+ g «on«5 + g g ^w^g, + 



[0100] Pre-multiplying (15) by (<f) T , where 1 S V £N and 
v*k, and using (1), (2), and (16) yields 



[0101] By (16), (18) and (19) we have determined 



[0094] Equating the coefficients of the 6G ; (i=l,2, . . . , m) 
terms in (14) yields 



[0095] Expanding z (1); k using normalized eigenvectors of 
(1) as basis vectors, one has 



4» = Z>W" 



[0102] Equating the coefficients of the SGfiGj terms in 
(14) yields 



2^ + ^4 J + ^4,= 

2'i \1 +\j +\, UZf„ +2'\f,,,U^ 

[0103] Expanding 
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[0104] using normalized eigenvectors of (1) as basis v< 
tors, one has 



4„ = Z^" 



-continued 



[0112] Expanding 

4* 



[0113] using normalized eigenvectors of (1) as basis vec- 
tors, one has 



[0106] is the coefficient of the u-th term and the number ir 
the parentheses in its subscript corresponds to that of 



4m = 2 rfvui-t" 
[0114] where 



[0107] Pre-multiplying (20) by ((j) k ) T and using (1), (2), 
(10), and (21) yields 



[0108] Substituting (21) into (11) and using (2) yields 



[0109] Pre-multiplying (20) by (f) T , where l^v^N and 
v*k, and using (1), (2), and (21) yields 



[0115] is the coefficient of the u-th term and the number ir 
the parentheses in its subscript corresponds to that of 



[0116] Pre-multiplying (25) by (<|> k ) T and using (1), (2), 
(10), and (26) yields 



|p4j ~A»i M Ai)ji -Ai)i M 42)« - A (i)i M 4)« 



[0110] By (21). (23) and (24) we have determined 



[0117] Substituting (26) into (12) and using (2) yields 



[0118] Pre-multiplying (25) by (f ) T , where l^v^N and 
v*k, and using (1), (2), and (26) yields 



[0111] Equating the coefficients of the 5G ; 5( ifiC I, t. 
(14) yields 



\dK dK dK 



3!(A> -A'/ ^dG^ 2 "' + dG~-^ m + acJ^W - A (i» Mz (2W " 
^4,' " 4^4* " 4^4. - A^.-M^y - A^M^,,] 

[0119] By (26), (28) and (29) we have determined 
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[0120] We proceed now to derive the perturbation solu- 
tions for the general p-th order terms in (5) and (6). Equating 
the coefficients of the SG^Gj . . . 5G s 5G t terms of order p in 
(14) yields 



[0126] The p-th order sensitivities of eigenvalues are 
obtained from (7) and (32). They depend on the eigenvalue 
and eigenvector sensitivities of orders up to p-2 and p-1 
respectively. Substituting (31) into (13) and using (2) yields 



= p !A*M4„ y ... „ + (p - m^Mjp-iv - s, 

A U)./ M 4-l)i ■ » + + A U)r M 4-l),j ■■■■ J + 
2 \(p - 2) ![Af 2w M4_ 21 ,.. a + A* 2)1 .,M4_ 2)J , 
A f:w M -i. IV,- ,] + ■■■ + P !A (p»j-- » M <** 



[0121] Expanding 



(4i, 1 ) 7 ' M 4-2u...» + • ■ ■ + (4w) TM 4-2)i/...J- ■ ■ + 

(p - 2) !2![(4_ 2 „ ,„ S ,) T M4 2)U + (4_ 2)J „ ) T M^ 2)jl + 

- + (4 P -2, i ,-...,) r M4 b >] + (p - D'[(4-i) J i...») 7 ' M 4«- + 

(4-d.i..J Tm 4) J + - + (4-iw.../ m 4)-]} 



[0127] Pre -multiplying (30) by (<)> V ) T , where 1 Sv^N and 
v*k, and using (1), (2), and (31) yields 



[0122] using normalized eigenvectors of (1) as basis v. 
tors, one has 



- + H«^i*...)-(p-D" 

4),j... » = X 'ow- <31) ( A (D. M 4-i)j...« +A (i)j M 4-i)i...» + ■■■ +A fi)r M 4-i; 

2 !(p - 2) !(A{ 2) , ( A"/4,, 2H...B + A{ :);/ 2) ■ + . . 



[0123] where 



A(,, i),...«,A"'4i).,' + ■■■ +A{,,-i, ! ,....,A':{i„) 



[0124] is the coefficient of the u-th term and the number ir 
the parenthesis in its subscript corresponds to that of 



[0128] By (31), (33) and (34) we have determined 



4* 



[0125] Pre-multiplying (30) by (c|) k ) T and using (1), (10), 
(31) and orthonorinalily relations of eigenvectors yields 




(32) 



^wMp-m- s) - 2 !(P " 2)!(A{ 2W M^ 2 „... „ + 

A (2),7 M 4-2)J - *, + ■■■ + A (2>» M 4-1W - 

(p-2)!2!(A^ 21i ... 1 ,M4 w + 

A t-2)>.- S ,M4).-| + - +A t-2W-.. r Mzf 2)B )] 



[0129] The p-th order sensitivities of eigenvectors can be 
subsequently obtained from (7). They depend on the eigen- 
value and eigenvector sensitivities of orders up to p-1. 

[0130] Equations (5) and (6) serve both the forward and 
inverse problems. In the former one determines the changes 
in the eigenparameters with changes in the stiffness param- 
eters. Damage detection is treated as an inverse problem, in 
which one identifies iteratively the changes in the stiffness 
parameters from a selected set of measured eigenparameters 
of the damaged structure: >. d k (k=l , 2, . . . , n,) and <p (l k (k= 1 , 
2, . . . , n^), where l = n x , n^N. Here X d k and <|) d k are 
assumed to be the perfect eigenparameters; simulated noise 
is included in the measured eigenparameters in the beam and 
frame examples that follow. Often we choose a set of n 
measured eigenparameter pairs to detect damage, i.e., 
%=n 4 ,=n. Let the number of the measured degrees of free- 
dom of (j) d k be N m ; N m =N and N m <N when we have 
complete and incomplete eigenvector measurements, 
respectively. With reduced measurements the unmeasured 
degrees of freedom of <|) d k is estimated first using a modified 
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eigenvector expansion method (see the beam and frame 
examples below) and <|) d k is mass-normalized subsequently. 
Only the component equations corresponding to the mea- 
sured degrees of freedom of <|) d k are used in (6). The system 
equations in Eqs. (5) and (6) involves n^n^f^ scalar 
equations with m unknowns, which are in general determi- 
nate if n x +n (|) N m =m, under-determined if n^+n (( ,N m <m, and 
over-determined if n x +n^N m >m. In the first iteration, X k and 
<|) k in (5) and (6) correspond to the eigenparameters of the 
structure with the initial stiffness parameters G ; ^ 0) . With the 
perturbation terms determined as shown above, the changes 
in the stiffness parameters 8G ; (1) , where the number in the 
superscript denotes the iteration number, can be solved from 
(5) and (6) using an optimization method discussed below. 
The estimated slii'i'ness parameters of the damaged structure 
are updated by G ; l 1) =G ; ,0) +8GA 1) . With the updated stiffness 
parameters, the eigenparameters, X k and <() k , in (5) and (6) are 
re-calculated from the eigenvalue problem (1) and the 
perturbation terms are re -evaluated. One subsequently finds 
6G ; (2) using the same method as that in the first iteration. 
This process is continued until the termination criterion, 
I'M 1; ! |<e, where L is the last iteration number and e is some 
small constant, is satisfied for all i=l,2, . . . , m. Note that 
after the w-th (l^w<L) iteration, we set G ; (w) to G hi if 
G; (w) >G hi , and to zero or some small stiffness value e G if 
Gj (w) <0. A flowchart for the iterative algorithm is shown in 
FIG. IB. When a single iteration is used, the iterative 
method becomes a direct method. 

[0131] FIG. IB shows steps included in a method for 
detecting structural damage in accordance with one embodi- 
ment of the present invention. The method includes as an 
initial step measuring one or more eigenparameters, k d k , <|> d k 
(Block 1). These eigenparameters are then compared with 
estimated eigenparameters associated with the stiffness 
parameters, Gi (0) (Block 2). The differences between the 
measured and estimated eigenparameters are then used by 
the perturbation method to establish system equations (5) 
and (6) (Block 3). The perturbation method may be a first or 
higher order multiple-parameter perturbation procedure. 

[0132] Next, an optimization method may be used to find 
the changes in the stiffness parameters 6G ; lw) (Block 4). 
These values are then compared to the predetermined value 
e (Block 5), and if the absolute values are less than e the 
stiffness parameters identified are set as G ; (w " l:i (Block 6). 
[0133] If the absolute values are not all less than e, the 
process proceeds along an iterative path where the stiffness 
parameters are first updated (Block 7). The stiffness param- 
eters are then bounded between 0 or e G and G h; (Block 8), 
and eigenparameters associated with the updated stiffness 
parameters are calculated (Block 9). The comparison of 
Block 2 is then performed based on these calculated, or 
estimated, eigenparameters. 

[0134] Optimization Methods 

[0135] Neglecting the residuals of order p+1 in (5) and (6) 
yields a system of polynomial equations of order p. When 
%+V^Sm. 8G ; (W) (i=l,2, . . . , m) at the w-th iteration can 
be solved from the resulting equations. There are generally 
an infinite number of solutions when n >t +n 4 ,N m <m, a unique 
solution when n x +n l|) N m =m and p=l, and a finite number of 
solutions when n x +n () ,N m =m and p>l. When n^+n^N^m, 
one generally cannot find 8G ; (w) to satisfy all the equations, 
and an optimization method can be used to find bG^ which 



minimize an objective function related to the errors in 
satisfying these equations. We use here the same notations, 
e x k and e^ k to denote the errors in satisfying the system 
equations in (5) and (6) respectively. Consider the objective 
function 



[0136] where W x k (k=l,2, . . . , n,) and W„, k (k=l, 2, . . . 
, n 0 ) are the weighting factors, and J is a function of 8G ; (w) 
when one substitutes the expressions for e x k and e^ in (5) 
and (6) into (35). When the first-order perturbations are 
retained in (5) and (6), the least-squares method can be used 
to determine 8G ; (w) which minimize (35) with W x k =W (| , k =l. 
Other weighting factors can be handled by dividing (5) and 
(6) by 




[0137] respectively. The method determines essentially 
the generalized inverse of the resulting system matrix, and 
is also known as the generalized inverse method. When the 
perturbations up to the p-th (p=l) order are included in (5) 
and (6), the gradient and quasi-Newton methods [25] can be 
used to determine 8G ; lw) iteratively. Unlike the generalized 
inverse method, the methods are applicable when other 
objective functions are defined. While the optimization 
methods are introduced for over-determined systems, they 
can be used when n^+n^N^m, in which case J=0 (i.e., 
e^e^O) when the optimal solutions are reached. 

[0138] Gradient Method 

[0139] To minimize the objective function in (35) at the 
w-th iteration, one can use the algorithm 



[0140] to update the changes in the stiffness parameters, 



[0141] a b >0 is the step size, and equals the gradient 




< dJ dJ dJ Y 
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[0142] associated with 



[0143] Note that the subscript b (b & 1) in all the variables 
in (36) denotes the number of nested iterations. The initial 
values used are 



[0144] The nested iteration is terminated when a b ||g b _ 
i|| 00 <Y, where is the infinity norm and y is some small 
constant, or the number of nested iterations exceeds an 
acceptable number, D. 

[0145] Quasi-Newton Methods 

[0146] Due to its successive linear approximations to the 
objective function, the gradient method can progress slowly 
when approaching a stationary point. The quasi-Newton 
methods can provide a remedy to the problem by using 
essentially quadratic approximations to the objective func- 
tion near the stationary point. The iteration scheme of these 
methods is given by (36) with f b .i=B b . 1 g b . 1 , where B^ is an 
approximation to the inverse of the Hessian matrix used at 
the b-th nested iteration, and the other variables the same as 
those previously discussed. Initially, we sel 



[0147] and B 0 =I, the identity matrix. The matrix Bt is 
updated using either the DFP formula 

(,)<;;;•„' - .,,>( - >r; , < 37 ) 
u^.'-ici;:!,, M a , -&,_!) 

[B b -i(g b -g b -i)][B b - l (g b -g b - l )] T 
(g b -g b -0 T B b -i(g b -g b -0 



the quasi-Newton methods, including the step size search 
procedure as outlined below, is shown in FIG. 2. 

[0150] Step Size Search Procedure 

[0151] The optimal step size for the gradient and quasi- 
Newton methods is determined in each nested iteration to 
minimize the function 



J(6&$_ l) -a b f b - l ) = F(a b ) 



[0152] with respect to a b . The search procedure is divided 
into two phases: an initial search to bracket the optimum a* b 
and a golden section search to locate a* b within the bracket, 
as shown in FIG. 2. 

[0153] Initial Bracketing. Choose the starting point Xj=0 
and an initial increment A>0. Let x^Xj+A, F 1 =F(x 1 ), and 
F 2 =F(x 2 ). Since for the gradient and quasi-Newton methods, 
f,,=g,, and it is along a descent direction of .1 when A is 
sufficiently small, one has F 2 <F-!. Rename 2A as A, and let 
x 3 =x 2 +A and F 3 =F(x 3 ). If F 3 >F 2 , stop and a* b is contained 
in the interval (x 1; x 3 ). Otherwise, rename x 2 as x 1 and x 3 as 
x 2 , then F 2 becomes F 1 and F 3 becomes F 2 . Rename 2A as A, 
and let x 3 =x 2 +A and F 3 =F(x 3 ). Compare F 3 and F 2 , and 
repeat the above procedure if F 3 <F 2 until F 3 >F 2 , with the 
final interval (x 1; x 3 ) containing a* b . 

[0154] Golden Section Search. If |x 3 -x 2 |>|x 2 -x 1 |, define a 
new point: 

x 4 =Ar 2+ 0.382(x 3 -x 2 ) (39) 

[0155] Otherwise, rename x 1 as x 3 and x 3 as x 1; and then 
define x 4 using (40). Let F 4 =F(x 4 ). If F 2 <F 4 , rename x 4 as x 3 , 
then I-'., becomes I-',. Otherwise, rename x, as x, and x 4 as x 2 , 
then F 2 becomes ¥ 1 and F 4 becomes F 2 . Compare |x 3 -x 2 | and 
|x 2 -Xj|, and repeat the above procedure until |x 3 -xj|f b _ 
ilU^ot' where e„ is some small constant. Then choose 




[0148] or the BFGS formula 

^ = ^ 1 + [l + ^-^f-^-^''J <38) 

I M - M , II M -'<■,• 

(6C$-6($. 1) )<g b -g b - 1 ?B b - l + 
B b - l (g b -g b - l )(g b -g b - l ) T 

[0149] The nested iteration is terminated when adlB^g,,. 
1 1 <y or the number of iterations exceeds D. A flowchart for 



MASS-SPRING SYSTEM EXAMPLE 

[0156] The algorithm discussed above is used to identify 
the stiffness parameters of a N-degree-of-freedom system 
consisting of a serial chain of masses and springs, such as the 
system shown in FIG. 3. Let the masses of the system be 
M ; =l kg (i=l, 2, . . . , N), and the stiffnesses of the 
undamaged springs be G h; =lN/m (i=l, 2, . . . , m), where 
m=N+l. The system is said to have a small, medium and 
large level of damage if the maximum reduction in the 
stiffnesses is within 30%, between 30 and 70%, and over 
70%, respectively. The mass matrix M is an NxN identity 
matrix, and the stiffness matrix K is a banded matrix with 
entries Ku-G^G^ (i=l,2, . . . , N), K^-K^— G 1+1 
(i=l, 2, . . . , N-l), and all other entries equal to zero. The 
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dK dK 
dG~i £md d~G^ 



[0157] have a unit value in entries (1, 1) and (N, N), 
respectively, and zero entries elsewhere. The nonzero entries 
of the matrices 




Since vanishing stiffness in any spring other than the two 
end ones in FIG. 1 can result in two decoupled subsystems, 
when G ; (w) <0 we set G ; (w) to e G =0.1N/m in the first two 
iterations and to O.lN/m in the remaining iterations for all 
the cases considered here. A relatively large value is 
assigned to e G in the initial iterations to avoid close eigen- 
values in the mass-spring system and small denominators in 
(19). This improves convergence especially when a small 
number of eigenparameter pairs are used. A smaller value is 
used for e G in the later iterations to improve the accuracy of 
stiffness estimation when there is a large level of stiffness 
reduction. Using the first-order perturbations and different 
numbers of eigenparameter pairs, the maximum errors in 
estimating the stiffnesses of the damaged system at the w-th 
iteration, defined by 



[0158] We look at a forward problem first with N=3 and 
m=4. The stiffnesses of the damaged system are G dl =G d3 = 
lN/m, G d2 =0.3N/m and G d4 =0N/m. The undamaged system 
is considered as the unperturbed system and the damaged 
system as the perturbed system. Based on the eigenparam- 
eters of the undamaged system, the eigensolutions of the 
damaged system are obtained using the first-, second-, and 
third-order perturbations, as shown in Table 1 below. The 
results show that even with the large changes in stiffness, the 
third-order perturbation solutions compare favorably with 
the exact solutions for the damaged system. The higher- 
order perturbation solutions can be used for large order 
systems when their direct eigensolutions become costly. 



\Gr-G d ,\ (40) 

= | nl ' L > Q~ d 

[0160] are shown in FIGS. 4A and 4B for all the itera- 
tions. In FIG. 4A p-1 and n- 1,2,3, and in FIG. 4B p=l and 
n=4,5, ... ,9. When n=l, the error decreases slow ly, though 
monotonically, and there is an estimation error of 1.5% at the 
end of iteration. While the errors can increase with the 
iteration number before approaching zero for n=3, they 
decrease monotonically for n=2 and n^4. All the stiffnesses 
are exactly identified at the end of iteration when n^2. Note 
that the number of the system equations equals and exceeds 
the number of unknowns when n=l and n^2, respectively. 
Since the system equations are linear, they have a unique 



TABLE 1 



1 0.165161 ("0.50000 1 ("0.285231 

0.65733 1 J 0.70711 I i 0.68836 I 

I 0.73528 J \ 0.50000 J [ 0.74129 J 

0.95541 1 f 0.707111 ( °- 95459 ] 

0.07840 I \ -0.00000 \ J 0.10607 1 

-IJ.2S47" J 1-0.71)711 J [-0.45962 J 



0.296351 (• 0.264451 
-0.74132 1 J -0.74875 1 
0.62482 J { 0.62362 J 



[0159] Consider now the damage detection problem with 
N=9, m=10, G d5 =0.5N/m, G d8 =0.7N/m, G dlo =0.8N/m and 
G di =lN/m (i=l,2,3,4,6,7,9). We set W^W^-l, e=0.001, 
Y=10" 10 , a ;-1 for all i, and D=500; the actual numbers of 
nested iterations in all the cases are much smaller than D. 



solution when n=l, and J has a unique minimum when n^2. 
With the small y the gradient method and the quasi-Newton 
methods using the DFP and BFGS formulas yield exactly the 
same results as the generalized inverse method (not shown 
here). Because the generalized inverse method does not 
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involve any nested iteration, it is the most efficient one 
among the four methods. While not shown here, the results 
indicated that that the quasi-Newton methods converge 
faster than the gradient method and the BFGS method has 
the similar performance to the DFP method. In what follows 
the BFGS method will be used with the higher-order per- 
turbations. With the second-order perturbations the errors 
shown in FIG. 4C decrease monotonically for all n (n=l,2, 
. . . , 9). The errors at w=l in the expanded view decrease in 
the order n=3,2,4,9,7,8, with the lines for n=5 and 7 virtually 
indistinguishable. In FIG. 4(A-C) the following symbols are 
used for different n: -—407 n=l; n=2; *-n=3; 

n=4; — n=5; — n=6; — n=7 - n=8; •• n=9. While 
the use of the second-order perturbations improves the 
accuracy of stiffness estimation in each iteration and reduces 
the number of iterations, it takes a much longer time to 
compute the higher-order perturbations and the associated 
optimal solutions. 

[0161] When only the first few eigenvalues are used, for 
instance, n x =5 and n^=0, the stiffnesses identified with the 
first-order perturbations, 



G< 4 >=(0.8" 



[0162] where the number in the superscript denotes the 
last iteration number, correspond to those of a different 
system with the same eigenvalues for the first five modes as 
the damaged system. The same stiffnesses are identified with 
the second-order perturbations. Similarly, when the first 
eigenvector is used, i.e., 11^=1 and n x =0, the stiffnesses 
identified with the first-order perturbations are those of a 
different system with the same eigenvector for the first mode 
as the damaged system: 



G< 9 >=(0.96 



(42) 

[0163] With the second-order perturbations the stiffnesses 
of the damaged system are identified. The stiffnesses iden- 
tified are not unique because the system equations in each 
iteration are under-determined. The solution given by the 
generalized inverse method here in each iteration is the 
minimum norm solution. Increasing the number of eigen- 
parameters used can avoid this problem. 

[0164] If the system has a large level of damage, i.e., 
G d5 =0.3 N/m, G dlo =0.1 N/m, with the other parameters 
unchanged, the stiffnesses of the damaged system are iden- 
tified with the first-order perturbations after 55 iterations 
when n=l and 6 iterations when n=2 and 3 as shown in FIG. 
5A. For n=4,5,6,7, the errors decrease monotonically and the 
number of iterations is reduced slightly, as shown in FIG. 
5B, which has an expanded view near w=l. With the 
second-order perturbations the errors shown in FIG. 5C 
decrease monotonically for n=l,2, . . . , 5, and the number 
of iteration for n=l is reduced from 55 in FIG. 5A to 4. The 
errors at w=l in the expanded view in FIG. 5C decrease in 
the order n=3,2,4,5, with the lines for n=2 and 4 virtually 
indistinguishable. The symbols used in FIG. 5A-C for n=l, 
2, . . . , 7 are the same as those in FIG. 4A-C. 

[0165] Finally, consider a large order system with a large 
level of damage: N=39, m=40, G dl2 =0.7N/m, G dl9 =G d37 = 
0.1 N/m, G d28 =0.8 N/m, G dl =l N/m (l=ii=i40 and i* 12,19, 
28,37), and the other parameters are the same as those in the 
previous example. With the first-order perturbations the 
e identified after 57 iterations when n=l, 



as shown FIG. 6A. Using a larger number of eigenparameter 
pairs (n=2,3,4,5) significantly reduces the number of itera- 
tions, as shown in FIG. 6B. The errors at w=2 in the 
expanded view in FIG. 6B decrease in the order n=4. 2.3,5. 
With the second-order perturbations the exact stiffnesses are 
identified after 6 iterations when n=l and 3 iterations when 
n=2, as shown in FIG. 6C, which has an expanded view for 
1S W £4. The symbols used in FIG. 6(A-C) for n=l, 2, . . . 
, 5 are the same as those in FIGS. 4A-C and 5A-C. 

FIXED-FIXED BEAM EXAMPLE 

[0166] The algorithm discussed above may be applied to 
detecting structural damage in an aluminum beam with fixed 
boundaries. The beam of length L t =0.7 m, width W=0.0254 
m, and thickness H=0.0031 m has an area moment of inertia 



— Wlf = 6.3058x 10 11 



[0167] and a mass density p=2715 kg/m 3 . The finite 
element method is used to model its transverse vibration. 
The beam is divided into N e elements, as shown in FIG. 7, 
with the length of each element being 



[0168] There are N e +1 nodes. With V ; and 9 ; denoting the 
translational and rotational displacements at node i (i=l, 2, 
. . . , N e +1), the displacement vector of the i-th (i=l,2, . . . 
, N e ) element is [Vi,e;,V ;+1 ,6 i+1 ] T . The Young's modulus is 
assumed to be constant over each beam element and that of 
the i-th element is denoted by G ; . The Young's modulus of 
the undamaged beam is G h =69xl0 9 N/m 2 . Hence G h; =G h for 
i=l,2, . . . , m, where m=N e . Small to large levels of damage 
correspond to reductions of the moduli tie lined for the 
Mass-Spring Example. The mass and stiffnes 
the i-th beam element are 



54 -13/,. 156 221, 



[0169] Using the standard assembly process yields the 
2(N e +l)x2(N e +l) global mass and stiffness matrices. Con- 
straining the translational and rotational displacements of 
the two nodes at the boundaries to zero yields the NxN M 
and K matrices, where N=2(N e -l) is the degrees of freedom 
of the system. The displacement vector of the system, 
involving the displacements of the 2 nd through N e -th node, 
is [V 2 ,e 2 ,V 3 ,e 3 , . . . ,V Ne , e N J T . The matrix 
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3G; 



[0170] (i=l, 2, . . . , m) can be obtained from K by setting 
G~l and G 1= . . . =G i _ 1 =G ;+1 = . . . =G N =0. The parameters 
W x k , W^, e, y, cr ; , and D are set to the same values as those 
previously discussed, and e G is set to 0.15 G h in the first two 
iterations and to 0.05 G h in the remaining iterations. The 
first-order perturbations are used unless indicated otherwise. 

[0171] Consider first the cases with N e =m=10 and N-18. 
When the system has a medium level of damage: 

G d =(l, 1, 1, 1, 0.5, 1,1, 0.7, 1, 0.8) T xG„ (44) 

[0172] the stiffness parameters of the damaged system are 
identified after 6 iterations with n=l. When the system has 
a large level of damage: 

G d =(l, 1, 1, 1, 0.3, 1, 1, 0.7, 1, 0.1) T xG h (45) 
[0173] the stiffness parameters of the damaged system are 
identified after 7 iterations with n=l. Consider next the cases 
with N e =20, 40 and 80. For the systems with medium and 
large levels of damage, the stiffness parameters of the first 10 
elements are given by (44) and (45), respectively, and those 
of the remaining elements are G h . In all the cases the 
stiffness parameters of the damaged systems are identified 
within 10 iterations when n=l. The numbers of iterations are 
reduced slightly when the second-order perturbations are 
used. Note that the system equations are over-determined 
when n=l. 

[0174] When only the translational degrees of freedom of 
an eigenvector are measured, a modified eigenvector expan- 
sion method is used to estimate the unmeasured rotational 
degrees of freedom. To this end, (() d k is partitioned in the 
form <|> d k =[(<|> dm k ) T , ('I'du 1 ')^' where (|) dm k and <|> dlJ k are the 
measured and unmeasured degrees of freedom of <|> d k , 
respectively. Similarly, <|> k in (6) is partitioned in the form 
( l )k =[( < l ) m k ) T ' ( ( l ) u k ) T ] T ' where (|) m k and (|> u k correspond to the 
measured and unmeasured components of (() d k , respectively. 
Since <|) dm k , <|) m k and (|) u k are known in each iteration, (|> du k is 
estimated from td u k= k[(<|>m k ) +( l ) dm k ] ( l ) u k > where the super- 
script+denotes generalized inverse. Once the rotational 
degrees of freedom of <|) d k are determined, (|> d k and <\> k are 
converted to their original forms and (() d k is mass-normal- 
ized. Only the component equations corresponding to the 
measured degrees of freedom of (|) d k are used in (6) and the 
system equations are determinate when n=l. The exact 
stiffness parameters of the damaged systems considered 
above can be identified. For the 10- and 20-element beams 
with the medium levels of damage, the stiffness parameters 
are identified after 6 and 24 iterations, respectively, when 
n x =l and n^=2. For the 10- and 20-element beams with the 
large levels of damage, the stiffness parameters are identified 
after 9 iterations when n x =l and n 4) =3 and 10 iterations when 
n=2, respectively. 

[0175] Finally, the effects of measurement noise on the 
performance of the algorithm are evaluated for the 10-ele- 
ment beam with the large level of damage. Simulated noise 
is included in the measured eigenparameters: 

k d k =k* d k +vR } k k* d k , WH^+v/VVd" (46) 



[0176] where X* d k and (|)* d k are the k-th perfect eigenvalue 
and eigenvector, respectively, R x k is a uniformly distributed 
random variable in the interval [-1,1], is a diagonal 
matrix whose diagonal entries are independently, uniformly 
distributed random variables in the interval [-1,1], and 
ve[0,l] is the noise level. Note that R x k and R 4) k are gener- 
ated for each measured mode. Each random parameter is 
generated 10 times and the average is used. Three different 
noise levels are considered: v=5'<, W'< and 20'r. When all 
the degrees of freedom of an eigenvector are measured, the 
stiffness parameters identified with n=l, 2, and 3 are shown 
in FIGS. 8A, 8B, and 8C, respectively. The following 
symbols are used for different noise levels: s ; v=5%; Q , 
v=10%; B, v=20%. When only the translational degrees of 
freedom of an eigenvector are measured, the eigenvector 
expansion method described above is used and the stiffness 
parameters identified with n=2 and n=3 are shown in FIGS. 
9A and 9B, respectively. The same symbols as those in FIG. 
8A-C are used in FIG. 9A-B for different noise levels. The 
stiffness parameters corresponding to v=0 in FIGS. 8 and 9 
are the exact values. It is seen that in the presence of noise, 
the stiffness parameters can be accurately identified with an 
increased number of measured eigenparameters. 

[0177] Space Frame Example 

[0178] The damage detection method can be applied to 
more complex structures, such as, for example, the modular, 
four bay space frame shown in FIG. 10. In this example and 
all the following examples the first order perturbation 
approach along with the generalized inverse method is used. 
The four-bay space frame example shown in FIG. 10 is 
made of extruded aluminum with 48 members, and is 8' tall, 
1.46' wide and 1.8' deep. The horizontal and diagonal 
members have the same cross-sectional dimensions (25.4 
mmxl2.7 mm) and the vertical members are "L"-angles 
with cross-sectional dimensions 50.8 mmx50.8 mmx6.35 
mm (thickness). All the members are connected through 
bolted joints. The frame is assumed to be fixed to the ground. 
The finite element method is used to model the 3-dimen- 
sional vibration of the frame, with each member modeled 
with four 12-degrees-of-freedom beam elements. The total 
degrees of freedom are 960 when the boundary conditions 
are applied. The Young's modulus is assumed to be constant 
over each member and that of the i-th (1 £is=m=48) member 
is denoted by G ; . The Young's modulus of the undamaged 
member is G h =69xl0 9 N/m 2 and G h; =G h . A vertical (mem- 
ber 1) and a diagonal (member 48) member in the first bay 
are assumed to have 70% and 30% reductions in Young's 
modulus, respectively, and all the other members are 
assumed to be undamaged. The first two vibration modes, 
corresponding to the bending of the frame along the x and 
y directions, respectively, are used to detect damage. The 
translational degrees of freedom of the 16 nodes, denoted by 
"S" in FIG. 10, along the x and y directions are assumed to 
be measured with 10% measurement noise. All the other 
degrees of freedom are estimated using the eigenvector 
expansion method discussed above and the measured modes 
are mass-normalized. The Young's moduli of all the mem- 
bers of the damaged truss are identified as shown in FIG. 10, 
with the maximum estimation error less than 7%. Note that 
the truss has closely spaced vibration modes. With the 
modes arranged in the order of increasing frequencies in 
each iteration, the method has effectively handled mode 
switching that has occurred in the damage detection process. 
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[0179] In summary, the damage detection method identi- 
fies stiffness parameters in structures, which have a small, 
medium, and large level of damage if the maximum reduc- 
tion in the stiffnesses is within 30%, between 30 and 70%, 
and over 70%, respectively: A large level of damage is 
studied in many examples because this poses the most 
challenging case, with sever mismatch between the eigen- 
parameters of the damaged and undamaged structures. 

[0180] The damage detection method as embodied and 
broadly described herein can be applied to structures that 
can be modeled with beam elements. A beam element is an 
element that has one dimension that is much longer than the 
other two. This element is very good at modeling 'T'-beams, 
rectangular beams, circular beams, "L"-angles, "C"-chan- 
nels, pipes, and beams with varying cross sections. Struc- 
tures that can be modeled with this element include, but are 
not limited to, lightning masts, light poles, traffic control 
poles, pillar type supports, bridges, pipelines, steel building 
frameworks, television, radio, and cellular towers, space 
structures, cranes, pipelines, railway tracks, and vehicle 
frames. Structures that can be modeled as beam elements are 
used simply for ease of discussion, and it is well understood 
that the damage detection method discussed above may be 
applied to structures that can be modeled by other elements 
using the finite element method or modeled by using other 
methods. 

[0181] Damage Detection Using Changes of Natural Fre- 
quencies: Simulation and Experimental Validation 

[0182] I 'or structures such as beams and lightning masts in 
electric substations, using only the changes in the natural 
frequencies can relatively accurately detect the localion(s) 
and extent of damage, even though the system equations are 
severely underdetermined in each iteration. This is an inter- 
esting finding as it is much easier to measure the natural 
frequencies than the mode shapes, and demonstrates the 
effectiveness of the iterative algorithm. Extensive numerical 
simulations on beams and lightning masts confirmed this 
finding. Experiments on the beam test specimens with 
different damage scenarios and a lightning mast in an 
electric substation validated the simulation results. The 
beam test specimens and the lightning mast are used as 
examples for demonstration purposes, and the method can 
be applied to other structures. Note that unlike the beam 
example shown earlier, where the Euler-Bernoulli beam 
finite element model is used, the Timoshenko beam finite 
element model is used in all the examples here. The Timosh- 
enko beam theory is found to be more accurate in predicting 
the natural frequencies of the lightning masts and circular 
beams than the Euler-Bernoulli beam theory. 

[0183] For a cantilever beam, simulation results show that 
the damage located at a position within 0-35% and 50-95% 
of the length of the beam from the cantilevered end can be 
easily detected with less than 5 measured natural frequen- 
cies, and the damage located at a position within 35-50% of 
the length of the beam from the cantilevered end and at a 
position within 5% of the length of the beam from the free 
end can be relatively accurately detected with 10-15 mea- 
sured natural frequencies. 



[0184] Numerical and Experimental Verification 

[0185] Cantilever Aluminum Beams 

[0186] Experimental damage detection results for four 
different scenarios are shown first, followed by various 
simulation results. 

[0187] Scenario 1: Evenly-Distributed Damage Machined 
from the Top and the Bottom Surfaces of the Beam Test 

[0188] The aluminum beam test specimen shown in FIG. 
12 is 45 cm long by 2.54 cm wide by 0.635 cm thick. It is 
divided into 40 elements (each element has a length of 1.125 
cm). The beam has a section (from approximately 10 cm to 
15 cm from the cantilevered end) of 5 cm long and 7.62E-4 
m thick machined both from the top and the bottom surfaces 
of the beam. This corresponds to 56% of damage (or 
reduction of bending stiffness EI) along the length of five 
elements (from the 9 th to the 13 rd element). Using the 
changes of the first 2 to 5 measured natural frequencies, 
damage is detected within 7 elements using 2 or 5 measured 
frequencies (from the 7 th to the 13 rd element with 5 mea- 
sured frequencies and from the 5 th to the 11 th element with 
2 measured frequencies) and within 8 elements using 3 or 4 
measured natural frequencies (both from the 6 th to the 13 rd 
element), as shown in FIG. 13. With 5 measured frequencies 
the average damage of the 7 elements in FIG. 13, from 6.75 
cm to 14.625 cm, is 46%. Note that the elements with 
damage less than 10% are not accounted for here. The extent 
of damage detected is slightly lower than the actual extent 
because the predicted damage occurs at 2 mote elements (the 
7 th and 8 th elements) than the actual one. The error results 
from the solution of the severely underdetermined system 
equations (5 equations with 80 unknowns). 

[0189] Scenario 2: The Same Aluminum Beam Test Speci- 
men as in Scenario 1 Clamped at the Other End. 

[0190] The same beam was tesled in the clamped-lree 
configuration with the clamped end reversed, which placed 
the damage from 25 cm to 30 cm (from the 23 rd to the 27 th 
element) from the cantilevered end. The damage detection 
results with 3 to 5 measured frequencies are shown in FIG. 
14. With 5 measured frequencies the average damage 
between 23.625 cm to 32.625 cm (from the 22 nd to the 29 th 
element) is 40%. Again the elements with damage less than 
10% are not included. 

[0191] Scenario 3: Undamaged Cantilever Aluminum 
Beam Test Specimen with the Same Dimensions as Above. 

[0192] An undamaged aluminum beam test specimen was 
clamped at one end with the same configuration as shown in 
FIG. 12. With the first 3 to 4 measured frequencies the 
maximum error for the estimated bending stiffnesses of all 
the elements is within 10%, as shown in FIG. 15. 

[0193] Scenario 4: A Cut of Small Width on a Cantilever 
Aluminum Beam Test Specimen, Shown in FIG. 16 with the 
Same Dimensions as Above. 

[0194] The beam shown below has a cut that is 0.4191 cm 
deep and 0.1016 cm wide, which corresponds to a 96'/? 
reduction in bending stiffness at the cut. The beam is divided 
into 45 elements and the cut is located in the middle of the 
23 rd element. With 3 to 5 measured natural frequencies, I he 
damage is detected at several elements surrounding the 23 rd 
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element, as shown in FIG. 17. It was found there was a 
roughly 50%, 40%, and 60% reduction in stiffness at the 
22'"', 23 rd , and 24 th element, respectively. The estimated 
damage extent at the above elements is lower than that at the 
cut because the damage is distributed along these elements. 
[0195] To examine the effectiveness and robustness of the 
damage detection algorithm, various simulations with dif- 
ferent damage scenarios were carried out. In this way, we 
can gain more insight concerning the accuracy of the finite 
element model, convergence of the estimated bending stiff- 
nesses of all the elements of the beam with the increased 
numbers of measured natural frequencies and/or mode 
shapes, and region of the beam within which the damage can 
be detected with few measured frequencies. With the can- 
tilever aluminum beam divided into 40 elements, the fol- 
lowing two simulations have the similar damage location 
and extent to those in Scenarios 1 and 2 in the experiments: 
[0196] Simulation 1: Uniform Damage Between 9.0 cm 
and 15.75 cm from the Cantilevered End of the Beam with 
the Same Dimensions as Discussed Above. 

[0197] Simulation 2: Uniform Damage Between 29.25 cm 
and 36 cm from the Cantilevered End of the Beam with the 
Same Dimensions as Discussed Above. 
[0198] Simulation results in FIG. 18 and FIG. 19 show 
that the damage can be relatively accurately detected with 
the first 3 to 5 measured frequencies, and can be accurately 
detected with the first 10 measured frequencies. 
[0199] With the cantilever aluminum beam divided into 
the same number of elements, two more simulations are 
presented here: one for a multiple damage scenario — 70% of 
damage at the 3 rd element and 30% of damage at the 20 th 
element, and the other for a 50% of uniform damage from 
the 16 th to the 18 th element (i.e., the damage is located at a 
position within 35-50% of the length of the beam from the 
cantilevered end). 

[0200] Simulation 3: Multiple Damage of the Beam with 
the Same Dimensions as Discussed Above: 70% of Damage 
at the 3 rd Element and 30% of Damage at the 20 lh Element. 

[0201] Simulation results in FIG. 20 show that the mul- 
tiple damage can lie relatively accurately detected with the 
first 3 to 5 measured frequencies, and accurately detected 
with the first 10 measured frequencies. 
[0202] Simulation 4: Uniform Damage from the 16 th to the 
18 th Element of the Beam with the Same Dimensions as 
Discussed Above. 

[0203] Simulation results in FIG. 21 show that the dam- 
age within 35-50% of the length of the beam from the 
cantilever end can be accurately detected with 10-15 mea- 
sured frequencies. 

[0204] Lightning Masts 

[0205] The lighting mast shown in FIG. 22 has two 

sections of constant cross-sections and an eccentric spike. 
The lengths of the lower and upper sections are 6.89 m and 
6.83 m, respectively, and that of the spike above the upper 
section is 1.4375 m. Both mode shapes and natural frequen- 
cies were measured for this mast. The mode shapes were 
measured by using a laser Doppler vibrometer. The natural 
frequencies were measured using both the laser Doppler 
vibrometer and an accelerometer. It takes about 0.5-1.5 days 



the mode shapes and about 30 minutes to 
measure the natural frequencies. It would be even more 
difficult to measure the mode shapes for some of the taller 
lightning masts, which are 100 and 130 feet tall. A finite 
element model was made using OpenFEM. Once the model 
was completed the measured and calculated natural frequen- 
cies and mode shapes were compared. The mast was 
expected to be undamaged, since by inspection the measured 
and calculated natural frequencies matched (Table 2). It was 
found that the first mode shape measurement is affected by 
wind and high modes are less affected by the wind. The 
measured frequencies are not affected by the wind and can 
be used for damage detection. 

TABLE 2 

Comparison of mcaMiicd and calculated (Irani the finite clement 



29.26 
45.32 
47.75 
54.10 
67.45 
74.70 



[0206] Damage detection was then performed using only 
the first 4 to 7 measured natural frequencies, as shown in 
Table 2. The mast is modeled by 40 elements; the lower 
section has 18 elements, the upper section has 15 elements, 
and the spike has 7 elements. The experimental damage 
detection results are shown in FIG. 23. The mast is modeled 
by 40 elements with 20 elements in each of the two sections. 
It appears that the elements around the joints near the middle 
of the mast (8.2 m in FIG. 23) and at the free end of the mast 
(13.72 m in FIG. 23) had some damage. This is most likely 
due to the fact that the joints were modeled in the finite 
element model as being infinitely stiff, which is almost never 
the case. Also, the spike is actually connected to the free end 
of the mast at two points. In the finite element model here 
it is assumed that the spike is attached to the mast along the 
whole line of contact that has a length of 0.34 m, which 
makes the model a little stiffer. Similar damage detection 
results were obtained with different numbers of measured 
frequencies. 

[0207] Simulations for the same lightning mast as shown 
in FIG. 22 with different damage scenarios are carried out. 
Shown in FIG. 24 are the simulation results for the mast 
with 70% damage between 0.76 m and 1.15 m from the 
ground. Shown in FIG. 25 are the simulation results for the 
mast with 50% damage between 6.89 m and 7.35 m from the 
ground. In both cases the location and extent of damage can 
be relatively accurately detected using the first 5 measured 
natural frequencies and accurately detected using the first 8 
or 10 measured frequencies. 
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[0208] Methods to Handle Some Ill-Conditioned System 
Equations 

[0209] When the first order perturbation approch is used, 
one needs to solve in each iteration a system of linear 
algebraic equations 

A6G=F (47) 

[0210] which are the linearized equations of (5) and (6), 
where A is the system matrix, F is the vector representing the 
differences between the measured and estimated natural 
frequencies and/or mode shapes, 8G is the optimal changes 
in the stiffness parameters to be found. While the gradient 
and quasi-Newton methods can be used, the generalized 
inverse method is most efficient because it does not involve 
nested iterations, and thus 6G=A + F, where A + is the gener- 
alized inverse of A. The problem (47) may be ill posed for 
certain cases, especially in the first several iterations, 
because A + can be relatively large and small changes in F 
can result in large changes in the solution 8G. Since the 
stiffness parameters cannot be negative or greater than the 
corresponding values of the undamaged structure, the solu- 
tion may not converge. 

[0211] Ill-conditioning problems do not occur in all the 
examples described above. Sometimes they can occur. Con- 
sider, for example, a cantilever aluminum beam of the same 
dimensions as those of the beam shown in FIG. 12. The 
beam is divided into 36 elements and assumed to have 30% 
of damage at the 5 th element and 90% of damage at the 18 th 
element. The translational degrees of freedom of the first or 
second mode shape vector are used to detect damage. The 
system equations are determinate, and the ill-conditioning 
problem occurs in the iterations. When a natural frequency 
is also used in additional to an incomplete eigenvector 
described above, the ill-conditioning problem can also 
occur. The following regularization methods can be used to 
handle the ill-conditioning problem: 

[0212] Method 1. Estimate A + from (A^+iiI)" 1 A T , 
where T| is a small positive constant that can be 
searched in several ways. One way is set T|=r|x 
1 .61 Sxniin( 10,||iSG|| ,). where , denotes the infin- 
ity norm, with the initial value T|=||A T A|| 00 , until 
IISGII^ is less than 0.8, for example. 

[0213] Method 2. To constrain the magnitude of 6G, 
we include it in the objective function to be mini- 
mized. For example, we can minimize the following 
objective function 



f(SG) = ]T f j (F l -A ij SGjf + f j (SG i f 



[0214] instead of (35), where N D =n x +n (|) N m is the number 
of equations in (5) and (6). The optimal solution can be 
obtained by using the generalized inverse method for the 
expanded system A*=[A;I] and the expanded vector F*= 
[F;0], where I is the mxm identity matrix and 0 is the mxl 

[0215] Note that the solutions from the regularization 
methods are not strictly the optimal solutions for the original 
objective function in (35). With accurate and sufficient 



measurement information and proper handling of the ill- 
conditioning problem, the system equations can become 
well conditioned in the last few iterations and regulation 
does not need to be applied. Consequently the stiffness 
parameters can be more accurately determined. Sometimes 
regulation may over constrain the magnitude of 6G, and the 
termination criterion ||8G|| co <e may be satisfied during the 
regulation process. In this case, we may set 



3.-0 </ ; 



[0216] to amplify the magnitude of 8G, and the iteration is 
terminated after the system equations become well condi- 
tioned. Since the second regulation method does not need to 
search for r\, it can be more effective when it works. When 
one carries out the damage detection procedures using 
different combinations of measured frequencies and mode 
shapes and obtains the same or similar stiffness parameters, 
one can have sufficient confidence on the results obtained. 

[0217] Using only the translational degrees of freedom of 
the first eigenvector and the first regulation method, the 
estimated bending stiffnesses of all the elements of the beam 
are shown in FIG. 26A. Using the translational degrees of 
freedom of the second eigenvector and the first regulation 
method, the estimated bending stiffnesses of all the elements 
of the beam are shown in FIG. 26B. In both cases the 
locations of damage are exactly detected and the extent is 
relatively accurately determined. Using the translational 
degrees of freedom of the first eigenvector along with the 
first natural frequency and the first regulation method, the 
bending stiffnesses of all the elements of the beam are 
exactly determined, as shown in FIG. 26C. 

[0218] Similarly, using the translational degrees of free- 
dom of the first eigenvector and the second regit lal ion 
method, the estimated bending stiffnesses of all the elements 
of the beam are shown in FIG. 27A. Using the translational 
degrees of freedom of the second eigenvector and the second 
regulation method, the estimated bending stiffnesses of all 
the elements of the beam are shown in FIG. 27B In both 
cases the locations of damage are exactly detected and the 
extent is relatively accurately determined. Using the trans- 
lational degrees of freedom of the first eigenvector along 
with the first natural frequency and the second regulation 
method, the bending stiffnesses of all the elements of the 
beam are exactly determined, as shown in FIG. 27C. 

[0219] Conclusions 

[0220] Thus, the sensitivities of eigenparameters of all 
orders may, for the first time, be derived using a multiple- 
parameter, general-order perturbation method. The higher- 
order solutions may be used to estimate the changes in the 
eigenparameters with large changes in the stiffness param- 
eters. The perturbation method may be combined with an 
optimization method to form a robust iterative damage 
detection algorithm. The gradient and quasi-Newton meth- 
ods can be used for the first or higher order system equa- 
tions, and the generalized inverse method can be used 
efficiently with the first order system equations because it 
does not involve nested iterations. Including the higher- 
order perturbations can significantly reduce the number of 
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iterations when there is a large level of damage. A modified 
eigenvector expansion method is used to estimate the 
unmeasured component of the measured mode shape. For 
many cases, the location(s) and extent of damage can be 
relatively accurately detected using only measured natural 
frequencies. Methods to handle ill-conditioned system equa- 
tions that may occasionally arise are developed and shown 
to be effective. Numerical simulations on different structures 
including spring-mass systems, beams, lightning masts, and 
frames show that with a small number of measured eigen- 
parameters, the stiffness parameters of the damaged system 
may be accurately identified in all the cases considered. 
Experiments on the different beam test specimens and the 
lightning mast in an electric substation validated the theo- 
retical predictions. The methodology can be readily applied 
to various operation structures of different sizes by incor- 
porating their finite element models or other mathematical 
models. 

[0221] One example of a practical application of this 
damage detection method is shown in FIG. 28, in which a 
structure 10 may represent any type of structure that can be 
modeled by beam elements, on which periodic damage 
assessment must be conducted. This includes, but is not 
limited to, structures such as lightning masts, utility poles, 
cell towers, and other such structures. A structure that can be 
modeled by beam elements is used simply for ease of 
discussion, and it is well understood that the general order 
perturbation method discussed above may be applied to a 
number of different structural members without departing 
from the spirit of the invention as embodied and broadly 
described herein. 

[0222] The structure 10 divided into elements E 1 through 
E n , and with nodes N ± through N n , is equipped with a sensor 
40, such as, for example, an accelerometer. The sensor 40 is 
configured to measure a response to an impact F or a force 
from a shaker applied to the structure 10. The impact F may 
be applied in the form of a single impact, as shown in FIG. 
28, or it may be in the form of a series of random impacts 
F-l through F n , as shown in FIG. 29. The location on the 
beam where the impact F is applied may be varied, but the 
impact F is preferably not applied at one of the nodal points 
of vibration modes whose natural frequencies and mode 
shapes need to be measured. In either instance, the sensor 40 
senses a response in the structure 10 to the impact F or the 
force from the shaker, and transmits that response to a 
processor 20 equipped with, among other elements, the 
damage detection method 30 discussed above. 

[0223] In accordance with the method as described above, 
upon receipt of the response signal from the sensor 40, the 
processor 20 accesses and performs the method 30 as 
previously discussed, and, as the data converges, determines 
an extent and location, by element, of any structural damage 
present in the structure 10. 

[0224] Although the convergence of the system to reso- 
lution is not dependent on where the impact F is applied to 
the structure 10, a signal to noise ratio of the response may 
be increased, and an efficiency of the system optimized, as 
the application point of the impact F is moved away from a 
fixed point N 0 , if one needs to measure the natural frequen- 
cies and/or mode shapes of lower modes. Similarly, although 
the system will converge regardless of the average energy 
and timing of the impact F or series of impacts F r F n applied 



to the structure 10, a signal to noise ratio of the response may 
be increased, and an efficiency of the system optimized 
based on a random series of impacts or the random shaker 
excitation to the structure 10, if the structure is large or 
slightly nonlinear or if there is an ambient excitation to the 
structure such as wind. 

[0225] The structure 10 may be excited in a number of 
different manners. For example, an impact F may be applied 
manually or a shaker may be used, at any given point on the 
structure 10, excluding the fixed point N 0 or the nodal points 
of vibration modes whose natural frequencies and/or mode 
shapes need to be measured. However, in an effort to 
increase the signal to noise ratio in the signal provided by the 
sensor 40 to the processor 30 if the structure is large or 
slightly nonlinear or if there is an ambient excitation to the 
structure such as wind and to provide the convenience and 
portability over the shaker test, a series of random impacts 
Fj-F n shown in FIG. 29 may be applied manually or by a 
specially designed device to generate the impact 1 . Although 
accurate and reliable damage assessments can be achieved 
regardless of how the impact F is applied to the structure 10, 
results may be improved using a random series impact 
method, which involves using a series of impacts F of 
random amplitudes and random arrival times. A series of 
random impacts has been shown to increase an energy input 
to the structure 10, improve the signal to noise ratio, 
especially in such situations as strong wind excitation, and 
average out slight nonlinearities that arise, for example, 
from bolted joints and extract linearized eigenparameters. 
[0226] FIG. 29 illustrates a structure 10, sensor 40 and 
processor 20 similar to those shown in FIG. 28. An impact 
F in FIG. 29 is applied to the structure 10, at a location other 
than N 0 or one of the nodal points of vibration modes whose 
natural frequencies and/or mode shapes need to be mea- 
sured, as a series of random impacts ¥ 1 through F n delivered 
at random amplitude and arrival time. These random impacts 
F^F,, may be delivered manually, or, as shown in FIG. 30, 
may he delivered in an automated fashion through the use of 
a random impact hammer device 1550. 
[0227] Different methods have been employed in conven- 
tional vibration testing in order to excite a test specimen. 
Shaker testing, in which a specimen is, simplistically, shaken 
in order to imparl a high level of energy, can produce a high 
signal to noise ratio, and can induce random excitation, 
which can average out slight nonlinearities and extract 
linearized eigenparameter parameters. However, shaker test- 
ing is not practically employed in the field on relatively large 
structures, and can be cost prohibitive to conduct. Single 
impact hammer testing addresses the shortfalls of shaker 
testing, in that it is portable and inexpensive to conduct. 
However, single impact hammer testing falls short where 
shaker testing is strong, in that low energy input of single 
impact hammer testing produces a low energy input, a low 
signal to noise ratio with no randomization. To address the 
need for a system which combines the advantages of shaker 
testing and single impact hammer testing, a Random Impact 
Series method for hammer testing is presented which yields 
a high energy, high signal to noise ratio, random system. A 
novel stochastic model is developed to simulate the random 
impact series produced manually and to generate a random 
impact series for a specially designed random impact device. 
[0228] FIG. 31 shows a random impact device 1550, 
according to one embodiment of the invention. Random 
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impact device 1550 includes random signal generating unit 
1560 and a random impact actuator 1570 with impact 
applicator 1580. The impact applicator 1580 preferably 
includes a sensor 1581, such as a force transducer, attached 
at its tip. Sensor 1581 is preferably configured to send data 
to the spectrum analyzer in stiffness parameter unit 103 in 
order to obtain mode shape information, as discussed above 
in connection with FIG. 1A. Sensor 1581 can be coupled to 
the spectrum analyzer in any manner know in the art 
including, but not limited to, wire, optical fiber or even a 
wireless connection. Random impact signal generating unit 
1560 generates outputs 1590 and 1591 to random impact 
actuator 1570. It should be appreciated that outputs 1590 and 
1591 could be coupled to the random impact actuator 1570 
in any manner known in the art including, but not limited, to 
cables, wires, optical fibers, wireless communications, and 
so forth. 

[0229] Output 1590 corresponds to the amplitude as 
discussed below. In particular, random signal generating unit 
1560 outputs a value corresponding to Uj ; , as discussed 
below. Output 1591 corresponds to x ; , as discussed below. 
[0230] Impact applicator 15S0 is shown, for purposes of 
illustration, as impacting structure 1600 in a reciprocating 
motion. Impact applicator 1580 is shown to have an impact 
path 1610. such that when a structure 1600 lies within the 
impact region 1610, impact applicator 1580 impacts struc- 
ture 1600 with a force of random amplitude that arrives at a 
random time, as discussed below. The impact applicator 
1580 and impact region 1610 as shown in FIG. 31 is one 
possible embodiment, and other shapes and impact regions 
may be used while still falling within the scope of the 

[0231] Although random signal generating unit 1560 is 
shown to output two outputs 1590 and 1591 that correspond 
to rj)j and Xj below, it should be understood that signal 
generating unit 1590 could output any signal or signals that 
ultimately results in random impact actuator 1570 driving 
impact applicator 1580 to impact structure 1600 with ran- 
dom arrival times x ; and random amplitudes r|i ; . 

[0232] As discussed above, vibration information is 
obtained by sensor 110, and is sent to stiffness parameter 
unit 103 via sensor coupler 1 13. The stiffness parameter unit 
103 sends stiffness parameters to the damage information 
processor 123. While the random impact device can be used 
for damage detection purposes as discussed above, it can 
certainly be used just for obtaining the natural frequencies 
and/or mode shapes of the structure for modal testing 
purposes. 

[0233] Experiments conducted on lightning masts by the 
applicant confirm that multiple impact testing performs 
better than single impact testing when there is wind excita- 
tion to the masts. Results are shown for a 65 foot tall mast 
with a 5 foot spike, as shown in FIG. 32, referred to 
hereinafter simply as a seventy foot tall mast. The mast has 
two constant cross section schedule 40 pipes of equal length. 
It can be seen from the results shown in FIGS. 33A-33B that 
the multiple impact test has much better coherence, is less 
noisy away from the resonances, and can pick up the modes 
that were missing in the single impact test. The multiple 
impact test is better at exciting the lower modes that can also 
be excited by the wind, which can be seen by comparing the 
frequency response functions (FRF) between 0 and 30 Hz. 



It can been seen from the FRF for the multiple impact tests 
that there are some modes that are missed at near 4 Hz and 
7 Hz with the single impact test, and the mode at 10 Hz is 
much improved compared to the single impact test. 

[0234] A stochastic model will now be discussed which 
describes the random impact series Fj-F^ modeled as a 
Poisson process. The Poisson process is one of a general 
class of processes that arise in problems concerning the 
counting of events in the course of time. The force pulses in 
the series are assumed to have an arbitrary, deterministic 
shape function and random amplitudes anil arrival times. 
The force signal in a finite time interval is shown to consist 
of wide sense stationary and non-stationary parts. The 
expectations of the average pow er densities associated with 
the entire force signal and the stationary part of the signal are 
derived and compared, and the power spectral density is 
related to the average power density and the autocorrelation 
function associated with the stationary part of the force 
signal. Numerical simulation is conducted to validate ana- 
lytical predictions, and a relationship between the Fourier 
transform and the discrete Fourier transform used in a 
numerical simulation is produced. Experiments on the four 
bay space frame, as shown in FIG. 10, validated the 
distributions of the random variables and the Poisson pro- 
cess. 

[1(235] Stochastic Model of a Random Impact Scries 

[0236] A random impact series is modeled here as a sum 
of force pulses with the same shape and random amplitudes 
and arrival times: 



[0237] where t is time, x(t) is the time function of the force 
signal, x ; e(0, t] is the random arrival time of the i-th pulse, 
and y(t-T[) is the deterministic shape function for all the 
pulses, is the random variable describing the amplitude of 
the i-th pulse, and N(t) is the number of the pulses that have 
arrived during the time interval (0,t] and is modeled as a 
Poisson process with stationary increments. All the pulses 
are assumed to be of width Ax, and y(t-X;) satisfies y(t-T,)=0 
if t<t; and t>x ; +Ax. 

[0238] Since a finite time record is used for the force 
signal in modal testing, we consider only the pulses that 
arrive during the time interval (0,T] of length T, as shown in 
FIG. 34. Note that y(0) and y(Ax) do not have to be equal 
to zero in this model. Because the pulses arriving after time 
t (t<T) will not affect the force signal at time t, the random 
process N(t) in (49) can be replaced by N(T). Also, if a pulse 
arrives at time T, it will vanish at time T+Ax. The last pulse 
in FIG. 34 arrives at time T N(T) and ends at a time between 
T and T+Ax. To include completely this last possible pulse, 
we consider the time interval t e(0, T+Ax]. Equation (49) is 
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mr> (50) . . mr, 

= g ^y(t - „), n e (0, r] and , e (0, T + Ar] J™ I £ ^ _ ^ I = ^ j£> + 



[0239] For the Poisson process N(t) with stationary incre- 
ments, the probability of the event {N(t)=n}, where n is an 



[0246] where F denotes Fourier transform, X(jw) is the 
Fourier transform of x(t), w is the angular frequency, and j = 
V-l. Let u=t-x ; , then t=u+x ; and dt=du. Equation (54) 



[0240] where X is the constant arrival rate of the pulses. By 
replacing t with T, the probability of the event {N(T)=n} is 



[0247] Since the integral in (55) is a deterministic func- 
tion, the expectation of the force spectrum is 



[0241] All the arrival times x ; , where i=l,2,3, . . . , N(T), 
are identically distributed, mutually independent random 
variables. Because the arrival rate of the pulses is constant, 
the arrival limes x,- arc uniformly distributed in [0. T] w ith 
the probability density function 



[0248] Using the expression for conditional expectations, 
E[E[V|U]]=E(V), where U=N(T) and 



[0242] Similarly, (i-1,2, . . . , N(T)) are identically 
distributed random variables, which are mutually indepen- 
dent and independent of the distribution of the arrival limes 
X;. While the distribution of r^ (i-1,2, . . . , N(T)) is not used 
in the subsequent derivation, it is assumed that i|> ; satisfy the 
Gaussian distribution with the probability density function 



[0249] we have from (56) 



E[XQoj)] = ^"^'j | WmjJ^Me-^i. 



[0243] where ,M=E[r[> ; ] is the mean and cr=L[i|' ,-]-E-[i|' ; ] 
is the variance of r[> ; . While the distribution in (53) is not 
used in the analytical derivations, it is used in the numerical 
simulation and validated experimentally for a manually 
applied random impacts on a four bay space frame, as shown 
in FIG. 10. 

[0244] Force Spectrum and its Expectation 

[0245] The force spectrum is the Fourier transform of the 

force signal in (50): 



[0250] Since Tj) ; e" iOT:i (i=l,2, . . . , n) are independent of 
each other and iJj ; is independent of e" JtOT ', we have 



[0251] Since (i=l,2, . . . , n) are identically distributed 
random variables and so are x ; (i=l,2, . . . , n), we have from 
(58) 
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[0252] Substituting (59) into (57) yields [0259] Let k=u-v, then u=k+v and du=dk. The double 

integral in (64) becomes 

E[XUco)] = < 6 °) 

£ P m (n, T)nE[^]^°e-^p T1 (tXt^t]^' yWe^du y( " )y(V)e " d " dV = 

[0253] Substituting (51) and (52) into (60) yields £ T dv jH^v + k)y(v)e-^dk 

ElXm] = Yj[ (XT> n< ]" £[ ^ l] [X (^"""H [0260] Interchanging the order of integration in (65) yields 

J ° J J ),(H)y(v) C -^»-"dHdv = 

£ T y(u)e-^d U j" T dk ^ T ~ k y{v ^k)y{v)e->^dv 

_ AE[»,](1 r y , u)e -^ du 
ju Jo 

[0261] Let v+k=v, then v=y-k and dy=dv . The first integral 
on the right-hand side of (66) becomes 

[0254] where Taylor expansion of e XT has been used. 
[0255] Average Power Density of x(t) in [0,T+Ax] and its 

Expectation ] dk i y(v+k)y(v)e-* J 'dv = 

[0256] The average power density of the force signal in r o r ^k 

(50) is denned as J J k J o y(7Mr-*)e" M <i7 

X(jai)X'(jco) (62) 
' S ' l(<,J) = — fTAr — [0262] Changing 7 in (67) back to v and substituting the 

resulting expression into (66) yields 

[0257] where X* (jw) is the complex conjugate of X(jw): 

j^' £ T y(u)y(y)e-**-* dudv = 

X'(jco) = g <p,e>^i f^ylule*™ du *'y(v + |*|)y(v)dve-^<*t 

[0258] Substituting (55) and (53) into (62) yields [° 263 ] Notin g that 



v|v + |A|).v(v)rfv 



[0264] is an even function of k, we have 

m y(v-i-\k\)y( v y-*" k dk = 



r 



y(y + \k\)y(v)co^k)dk 
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[0265] Substituting (69) into (68) and substituting the 
resulting expression into (64) yields 



4- |t|)jr(y)<fvCOKufatt 



-continued 

= XTE[^] + 2X 2 E-[^ i ] — 

[0271] Substituting (74) into (71) yields 



[0266] Since the double integral in (70) is deterministic, 
we have 



[0267] Since 

CI 



^,>,-sin|.jMr,--r,„)| =(1.<4 ; I/=1.2 V(7')) 



[0268] arc identicallv distributed random variables, and si 
are x ; (i-1, 2, . . . , N(T)), we have 



nv(7W) 1 nvm/vm 

£ Z ^ + Z Z - t») 



l/V(r )l A, : + [/v : (D-/v(r)] 1 A, 2 co- 



[0269] Since Xj and x 2 are independently, uniformly dis- 
tributed random variables in [0,T], the probability density 
function of X-j-x 2 is 



£[SiMl = TTI^f^f^ "" y{v + WMvJdvcost^d*] x (?5) 

{2A^[,/„ \ l "^'" n +.177:1*11} 

[0272] Mean and Autocorrelation Functions of x(t) 

[0273] The first-order cumulant function of x(t) in (50), 
K^xO)], is equal to its mean function, E[x(t)]. The second- 
order cumulant function of x(t), K 2 [x(tj)x(t 2 )], is related to 
its autocorrelation function, E[x(tj)x(t 2 )] through 

ElxihXQhKlxitM'^Mk^MQ] (76) 
[0274] where t 1 and t ? are any two time instants in [0,T+ 
Ax]. Following the derivations, the first- and second-order 
cumulant functions of x(t) are 



k- , |a(/) ]=/•;[ v(r) hX£[il>JJ 0 T y(t-a)da 



(77) 



[0275] where te[0,T+Ax]. Let t-a-u in (75), then da=-du. 
We have from (77) 



mn\ = \Wx\\ v(«V/« 



[0276] Let 

W(<HoM»>*< (80) 
[0277] Noting that y(u)=0 when u<0 and u>Ax, we have 



from (79) 



I T -A\t\ _ (73) (mtiWit), 0<f<Ar 

- r2 (r) = | ' T - T - T E[x(t)] = hE» l ]W(AT), At < t < T 

(,0, elsewhere ( A£[^,][W(Ar) - W(t - T)], T<t<T + . 



[0270] Using (51) and (73) in (72) yields 
/| v v , , | = y P „ T „ E[I : ] + 

^P |A ,|(n,r)(n 2 -n)£[^] 



[0278] where T>Ax is assumed. Let tj-o^u and t 2 -t 1 =k in 
(78), then da — du and t 2 -a=u+k. We have from (78) 



'il f l y(u)y(u- 
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[0279] When |k|>Ax, y(u)y(u+k)=0 and hence t 2 )= 
0. When QSk^Ax, we have from (82) for different t x 



[0284] Fourier Transform of the Mean Function of x(t) and 
its Equivalence to E[X(jw)] 

[0285] Applying the Fourier transform to E[x(t)] in Eq. 
(81) yields 



AE[iA? ]£ l y(u)y(u + k)du, 0 < h < At - * 
XE W\] k y(u)y(u + k)du, At - k < r, < T 



■k)du, T<h <T + £*.T-k 
T + At — k<t\ <T + 



[0280] When -Ax£k<0, we have from (82) for different t 2 



\E[if/\ ] P y(u)y(u + k)du, -k < h < At 

XEW\ ] f%My(" + k)du, At < ft < T - k 

\E[<P\ ] J"' y(u)y(u + k)du, T-k< h <T + 

o, o s n < -* 



[0281] Let u+k=v in (84), then u=v-k and du=dv. We have 
from (84) after changing v back to u 



F{E[x(t)]} = *Wi\{£ T W)e-**dt + f\(&T)e-*«dt + <88) 
£ +&T [W(M-W(t-T)]e-'"d^ 



[0286] The three integrals in (88) are referred to as 
and I 3 , respectively. Consider first the third integral ir 



h\ [W(AT)-W(t-T)]e-^'dt 



[0287] Let t-T=8 in (88), then t=T+6 and dt=d6. We have 
from Eq. (89) 



[0288] Changing 6 in (90) back to t and combining l ± and 
I 3 yields 



1 iT W(AT)e-'""<&+(l-e-' <u ' I )J 0 AT W(f)e- 



[0289] Adding I, to (91), simplifying the expression, and 
substituting it into (88) yields 



XE[^]p* k y(u-k)y(u)du, 
\E[^]£ T ^y( U -k)y( U )du, 
XE[^]j^y(, 



y(u-k)y(u)du, T — 



F{E[x(t)]} = 

m^-e-^fy^dr^jy^ 

e-*"dt + f T W{fr)e-*-*dtj = \E[i/ ll ] 

a-e-^i± + £\^-^d t ] 



[0282] Combining the second equations in (83) and (85), 
we have for t x and t 2 in [Ax,T] 



K xx ( h , h )=XE[^ l \^ T m y(u)y(u + 



[0283] Since by the second equation in Eq. (81), E[x(t)] is 
a constant for t e[Ax, T], and by (86), K xx (t-L, t 2 ) is a function 
of k=t 2 -t 1 for tj and t 2 in [Ax, T], x(t) is a wide-sense 
stationary random process in [Ax, T]. Substituting the sec- 
ond equation in (81) and (86) into (76) yields the autocor- 
relation function for t t and t 2 in [Ax,T] 



R xx (k)nE[x(hMh)] = 

XEJtii] K y{u)y(u + \k\)du+X 2 E l VI,iW 2 {bT) 



[0290] We will show that F{E[x(t)]] in (92) is equivalent 
to E[X(jw)] in (91). By (80), we have 



(-oo, 0) 
I(t)dt, re[0,AT] 



[0291] The integral in (61) can be v, 



f% {u) e^du = f\ (u) du + £ T y{u){ ^- m u 

= £ T y M du-J.£y W £e^d V d 
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[0292] Interchanging the order of integration in the double [0301] where X s (jc>)= | Xx 1 x(t)e |OJt dt. Taking the expec- 
integral in (94) and using (93) yields tation of (96) yields 



jTW*-*.= (95) ^ = ¥^£1 

= W(A,){ 1 + ^ 4r [|^-l>-^v} 



[0302] Let k=t 2 -t 2 in (97), then dk=dt 2 . Since 
E[x(t-!)x(t 2 )]=R xx (k), (97) becomes after interchanging the 
order of integration 



[0293] Substituting (95) into (61) yields (91). This shows 
that the expectation E and the Fourier transform F are 
commutative as both are linear operators. 

[0294] Equation (92) consists of two parts: the first part, 



AE[lAl](l-«e-"'">Vmr)^ 



[0303] Substituting (87) into (98) yields 
E[S 2 M]=AEWfi 



[0295] is the Fourier transform of the stationary part of 

E[x(t)] in [Ax,T], and the second part, 2A 2 £ 2 [iA,]w 2 (Ar)J^ " ' (l - j^) eMd k 

A£[^](l - £ ->"<)W(Ar) - l) c -**dr, [0304] Noting that y(u)y(u+|k|)=0 when |k|>Ax and 

[0296] is the sum of the Fourier transforms of the nonsta- X yb + W>yW>^ - 

tionary parts of E[x(t)] in [0,Ax] and [T,T+Ax]. When Ax^O, 

since 

[0305] an even function of |k|, we have from (99) 



[0297] is finite, 



E[S 2 (u)} = 2A 2 E 2 Wi}W-y\j) Ar) + AZT[^f ] 



[0306] where T>2Ax has been assumed. 

[0307] The power spectral density of x(t) can be obtained 

from (100) by increasing T to infinity 



S x (a>) = KmS 2 (o>) = 7xX 1 E 1 [ip i ]W 2 {b.T)6(u) + ( m ) 



[0298] and consequently the second part of E[X(jw))], 
approaches zero. 

[0299] Average Power Density of x(t) in [Ax,T], Its Expec- A£[lA ' ] X/ %« + l*l)/(«)cosM:)<tt 

tation, and Power Spectral Density 

[0300] Since x(t) is stationary in [Ax,T], the average power [ 0308 ] where we have used 
density of x(t) in [Ax,T] is defined as 
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[0309] in which 8(0) is the Dirac delta function. The 
power spectral density can also be obtained from (98) by 
increasing T to infinity 

(103) 

[0310] where R xx (-k)=R xx (k) has been used. Equation 
( 103) is the well-known Wiener-Khintchine theorem, which 
states the power spectral density is the Fourier transform of 
the autocorrelation function. Substituting (87) into (83), and 
noting that I(u)I(u+|k|)=0 when |k|>Ar and the Fourier 
transform of 1 is 2;to((»), yields (101). Note that the power 
spectral density is only defined for a wide -sense stationary 
process with an infinite time record. When the mean ampli- 
tude of each pulse EfipJ is not equal to zero, there is an 
associated delta function in the power spectral density. 

[0311] Comparison of EfS^co))] and E[S 2 (w)] 

[0312] By (100), E[S 2 (w)] consists of two parts: the first 
part, 



[0317] When y(v)=W 0 5(v), where W 0 is : 
have from (104) 



y!v + |/:|iy!VK/icos(„/u/A = 



4 liin|j^ A V 0 5(v)£--»" v clv| = W 0 2 



[0318] Since |k|<Ax in (100), we have |k|-»0 as Ax^O and 
hence 



[0313] which depends on the arrival rate X, the mean 
amplitude of each pulse F[i|',], the total area of the normal- 
ized shape function W(Ax), and the time length T-Ax, 
describes the average effects of the stationary part of x(t) in 
[Ax,T] and is referred to here as the first-order statistical 
power density, and the second part, 



[0319] Substituting (105) and (106) into (100) and noting 
that W(Ax)=W 0 when y(v)=W 0 5(v), and substituting (105) 
into (75) and noting that T+Ax-»T as Ax-»0, yields 



,vUn,/„(l- 7 - I -jco^. 



£[5 1 M] = 2A 2 £ 2 [iA I ]- 



[0314] which depends on the mean square amplitude of 
each pulse EPi^ 2 ] and the shape function y(0) in addition to 
~k, T, and Ax, describes the variational effects of the station- 
ary part of x(t) and is referred to here as the second-order 
statistical power density. While EfS^cu)] in (75) consists of 
two parts, the first part, 



v(v + W,v<v,rfveo«^ «Ik - 



[0315] depends also on the shape function y(D) as the 
shape function is used in calculating the power associated 
with the nonstationary parts of x(t) in [0,Ax] and [T, T+Ax]. 

[0316] When the shape function is a delta function (i.e., 
Ax-»0), the nonstationary parts of x(t) vanish and E[S,((n)] 
in (75) can be shown to be equivalent to E[S,(c>)] in (100). 
By (68) and (69), we have 



[0320] By (104) the power spectral density in (101) can be 
expressed as 

S x («)=2 3t X, 2 iT 2 [^ 1 ]W 2 (AT)8((0) + A£[ 1 |. 1 2 ]|F(,-»)| 2 (108) 

[0321] where Y(j ro)= J 0 AT y(t)e- j<ot . When T^°o in (75), we 
have by using (104) 



.V(...i-/h» 7Wa./.|.s,r«„|-2.-T/.7. 1 v, l|) './""| V«..+>y. 
[1>iW)| 2 



(109) 



[0322] Equation (109) has a slightly different form from 
that of (108) because the power associated with the nonsta- 
tionary parts of x(t) is included in (109). When y(t)=W fJ 6(t) 
(108) and (109) reduce to 



(110) 



[0323] Equation (110) can also be obtained from (106) by 
letting T-»oo and using 
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EXAMPLES AND NUMERICAL SIMULATION 



[0324] When the shape function of the pulses is repre- 
sented by a half sine wave, i.e., 



shown as a solid line in FIG. 38, is calculated using (75) . 
The curve for 10 logjEfS^ffi)]} with X=l/s and the other 
parameters unchanged is shown as a dashed line in FIG. 38. 
It is seen that EfS^co)] increases by 4.14 to 15.6247 times 
in the frequency range shown when X is increased from 
X^l/s to >t 1 =4.14/s. This result can be shown by using (75) 



[0325] where H(0) is the Heaviside function, we obtain by 
using (60), (64), and (103) 



E[X { jco)] = 



(oj 2 £iT 2 -n 2 ) 2 (T + . 



-2 A - /_-.<■> i »o- ./-I. I 



E[S 2 (co)] = AEt^Ar 2 - 



(w 2 Ar 2 -7T 2 ) 3 (7--Ai-) 



[+(8sin(wAr)w^ + 
27-cos(wAt)wV +2Tn 2 co 2 )Ar 2 ] 



fo 2 Ar 2 -^) 3 (T-Ar) 



Ajf^iAjr 2 +A l TE[if/ 2 l ] _ ( 
A;/.:|.',,|/:+.l : //.l*i| " 

A^^]r + A l£[ ^] 
A|£2[^]r + A 2 £[^] 



Hm£[S,(; W )]| A=Al 



[0327] This shows that a larger arrival rate X would 
increase the energy input to the structure over the entire 
frequency domain. 

[0328] Numerical simulation is undertaken next to vali- 
date the analytical predictions. The random number N(T) 
satisfying the Poisson distribution in (5 I ) with >.=4. 14 s and 
T=8 s is generated using MAT! AH. Similarly, the random 
numbers corresponding to the random variables %■ (i=l, 2, . 
. . , N(T)), satisfying the uniform distribution in (52), and the 
random numbers corresponding to the random variables t|) ; 
(i=l, 2, . . . , N(T)), satisfying the Gaussian distribution in 
(53) with ,m=0.8239 N and a 2 =0.7163-0.82392 N 2 =0.0375 
N 2 , are generated. Using the shape function constructed 
earlier, a sample function of x(t) in (50) at time t=rh, where 
r=0, 1, , 



)Ar : ll-cosfc>(7--AT)) 
n 2 oj 2 (T - At) 



[0326] Consider next the normalized shape function y(t) 
shown in FIG. 35 with unit maximum amplitude. It is 
obtained by averaging a series of normalized force pulses 
from impact tests on the four-bay space frame as shown in 
FIG. 10. There are 21 sample points in the shape function, 
which are connected, as shown in FIG. 34. Other parameters 
used are T=8 s, Ax=20xT/l 024=0. 15625 s where h=T/1024= 
0.0078125/s is the sampling interval, X=4.14/s, tf$i\= 
0.8239N, and EP^^OJieSN 2 . The curve for E[x(t)] in the 
time interval from 0 to 8.15625 s, shown as a solid line in 
FIG. 36, is calculated using (81). It is seen that E[x(t)] 
increases from 0 lo 0.0352 N in the first 0.15625 s, remains 
at 0.0352 N from 0.15625 s to 8 s, and decreases from 
0.0352 N to 0 in the last 0.15625 s. The curve for 20 
log|E[X(jo>)] in the frequency range from 0 to 50 Hz, shown 
as a dotted line in FIG. 37, is calculated using (61). The 
curve for 10 loglEfS^ai))]} in the same frequency range, 



[0329] denoted by x r , can be obtained. The discrete Fou- 
rier transform (DFT) of the time series {x r } is calculated 
using MATLAB. The DFT of the series {x r } is defined by 



".-it"-* 



[0330] where q=0, 1, . . . , R-l. Equation (118) is an 
approximate formula for calculating the coefficients of the 
Fourier series of a periodic function whose values in the 
period [0, T+Ax] are given by those of x(t): 
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[0331] The Fourier components X q correspond to harmon- 
ics of frequency 



[0332] Recall that the Fourier transform of x(t) in (50) is 
given by 



[0334] in (120) and compare the resulting expression with 
(119), we find that X q in (118) multiplied by T+Ax provides 
an approximate value of X(jco) at frequency co q . Similarly, 
X X* q multiplied by T+Ax provides an approximate value of 



[0335] at frequency co q . By averaging 5000 sample series 
of {x r }, we obtain the curve for E[x(t)], shown as a dotted 
line in FIG. 36. By averaging 1000 sample series of {20 
log[( T+At)|XJ]], we obtain the curve for 20 log|li[X(jc.))]|, 
shown as a solid line in FIG. 37. By averaging 100 sample 
series of {10 log[(T+Ax)|X q X* q |]}, we obtain the curve for 
10 logjEfS^cu)]}, shown as a dotted line in FIG. 38. The 
numerical results are in good agreement with the analytical 



[0336] The stochastic model was experimentally validated 
for an experimenter conducting manually a random series of 
impacts on the four bay space frame as shown in FIG. 10. 
One hundred ensemble averages were used. The experimen- 
tal probability density functions of the arrival time, the 
number of arrived pulses, and the pulse amplitudes are in 
good agreement with the analytical values, as shown in 
FIGS. 39-41. 

[0337] Thus, the system and method for detecting struc- 
tural damage and the random impact series method as 
embodied and broadly described herein can be applied to an 
unlimited number and type of structures to provide auto- 
mated, reliable damage detection and assessment and to 
conduct modal testing. This system could be further auto- 
mated to conduct periodic tests and provide results to a 
centralized monitoring section. Regular health monitoring of 
these types of structures could provide additional protection 
against potential failure, as well as a characterization of 
usage and wear over time in particular environmental con- 
ditions for predicting useful service life. 
[0338] The foregoing embodiments and advantages are 
merely exemplary and are not to be construed as limiting the 



present invention. The present teaching can be readily 
applied to other types of systems. The description of the 
present invention is intended to be illustrative, and not to 
limit the scope of the claims. Many alternatives, modifica- 
tions, and variations will be apparent to those skilled in the 
art. In the claims, means-plus-function clauses are intended 
to cover the structures described herein as performing the 
recited function and not only structural equivalents but also 
equivalent structures. 
What is claimed is: 

1. A system for determining stiffness parameters of a 
structure, comprising: 

a sensor arranged to measure vibrations of said structure 
and output vibration information; and 

a stiffness parameter unit for receiving said vibration 
information, determining natural frequency data of said 
structure, and determining the stiffness parameters of 
said structure using said natural frequency data. 

2. The system according to claim 1, further comprising 
multiple sensors arranged to measure vibrations of said 
structure and output vibration information. 

3. The system according to claim 1, wherein said stiffness 
parameter unit comprises an iterative processing unit. 

4. The system according to claim 1, wherein said stiffness 
parameter unit comprises an outer iterative processing unit 
and an inner iterative processing unit. 

5. The system according to claim 3, wherein said iterative 
processing unit determines said stiffness parameters using a 
first order perturbation process. 

6. The system according to claim 3, wherein said iterative 
processing unit determines said stiffness parameters using a 
higher order perturbation process. 

7. A system for determining stiffness parameters of a 
structure, comprising: 

a sensor arranged to measure vibrations of said structure 

and output vibration information; and 
a stiffness parameter unit for receiving said vibration 

information and determining said stiffness parameters 

with an iterative processing unit. 

8. The system according to claim 7, wherein said iterative 
processing unit comprises an outer iterative processing unit 
and an inner iterative processing unit. 

9. The system according to claim 7, wherein said iterative 
processing unit determines said stiffness parameters using a 
first order perturbation process. 

10. The system according to claim 7, wherein said itera- 
tive processing unit determines said stiffness parameters 
using a higher order perturbation process. 

11. A stiffness parameter unit for determining stiffness 
parameters for a structure, comprising: 

an input for receiving vibration data related to the struc- 



i processing unit for receiving said spectral 
data and outputting said stiffness parameters using 
natural frequencies of the structure. 
12. A stiffness parameter unit for determining stiffness 
parameters for a structure, comprising: 

an input for receiving vibration data related to the struc- 
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an analyzer for converting said vibration data to spectral 
data; and 

an interative processing unit for receiving said spectral 
data and outputting said stiffness parameters using a 
perturbation process. 

13. The stiffness parameter unit according to claim 12, 
wherein said perturbation process comprises a first order 
perturbation process. 

14. The stiffness parameter unit according to claim 12, 
wherein said perturbation process comprises a higher order 
perturbation process. 

15. A system for determining damage information of a 
structure, comprising: 

a sensor arranged to measure vibrations of said structure 
and output vibration information; 

a stiffness parameter unit for receiving said vibration 
information, determining natural frequency data of said 
structure, and determining the stiffness parameters of 
said structure using said natural frequency data; and 

a damage information processor for receiving said stiff- 
ness parameters and outputting damage information. 

16. The system according to claim 15, wherein said 
damage information processor outputs damage location 
information or extent of damage information. 

17. A system, comprising: 
a structure; 

a sensor arranged to measure vibrations of said structure 
and output vibration information; and 

a stiffness parameter unit for receiving said vibration 
information, determining natural frequency data of said 
structure, and determining the stiffness parameters of 
said structure using said natural frequency data. 

18. The system according to claim 17, further comprising 
a damage information processor for receiving said stiffness 
parameters and outputting location of damage. 

19. The system according to claim 18, wherein said 
damage information processor comprises a damage location 
processor for determining damage location information. 

20. The system according to claim 18, wherein said 
damage information processor comprises a damage extent 
processor for determining extent of damage information. 

21. The system according to claim 18, wherein said 
damage information processor comprises a damage extent 
processor for determining extent of damage information and 
a damage location processor for determining damage loca- 
tion information. 

22. The system according to claim 17, wherein said sensor 
comprises a velocimeter. 

23. The system according to claim 17, wherein said sensor 
is attached to said structure. 

24. The system according to claim 17, wherein said sensor 
is not attached to said structure. 

25. The system according to claim 17, wherein said 
stiffness parameter unit further comprises a spectral ana- 

26. The system according to claim 17, wherein said 
structure comprises a beam. 

27. The system according to claim 17, wherein said 
structure comprises a truss. 

28. The system according to claim 17, wherein said 
structure has a longest dimension less than 1.5 meters. 



29. The system according to claim 17, wherein said 
structure has a longest dimension less than 2.5 meters. 

30. The system according to claim 17, wherein said 
structure has a longest dimension less than 10 meters. 

31. The system according to claim 17, wherein said 
structure has a longest dimension less than 50 meters. 

32. A device, comprising: 

a random signal generating unit for generating first and 
second outputs; 

a random impact actuator for receiving said first and 
second outputs; and 

an impact applicator coupled to said random impact 
actuator and having an impact region; 

wherein said random impact actuator drives said impact 
applicator such that the force and arrival times of said 
impact applicator at said impact region are random. 

33. The device of claim 32, wherein said random impact 
actuator drives said impact applicator in accordance with 
said first and second outputs. 

34. The device of claim 33, w herein the first and second 
outputs comprise independent random variables. 

35. The device of claim 34, wherein the first and second 
outputs determine the force and arrival times, respectively, 
of the impact applicator at the impact region. 

36. A system, comprising: 

a structure; 

a random impact device for inducing vibrations in said 
structure; 

a sensor arranged to measure vibrations of said structure 
and output vibration information; and 

a stiffness parameter unit for receiving said vibration 
information, determining natural frequency data of said 
structure, and determining the stiffness parameters of 
said structure using said natural frequency data. 

37. The system of claim 36, wherein the random impact 
device comprises: 

a random signal generating unit for generating first and 
second outputs; 

a random impact actuator for receiving said first and 
second outputs; and 

an impact applicator coupled to said random impact 
actuator and having an impact region; 

wherein said random impact actuator drives said impact 
applicator such that the force and arrival times of said 
impact applicator at said impact region are random. 

38. The device of claim 37, wherein said random impact 
actuator drives said impact applicator in accordance with 
said first and second outputs. 

39. The device of claim 38, wherein the first and second 
outputs comprise independent random variables. 

40. The device of claim 39, wherein the first and second 
outputs determine the force and arrival times, respectively, 
of the impact applicator at the impact region. 

41. A system for determining stiffness parameters of a 
structure, comprising: 
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a sensor arranged to measure vibrations of said structure 
and output vibration information; and 

a stiffness parameter unit for receiving said vibration 
information, determining mode shape information, and 
determining the stiffness parameters of said structure 
using said mode shape information. 

42. The system according to claim 41, further comprising 
multiple sensors arranged to measure vibrations of said 
structure and output vibration information. 

43. The system according to claim 41. wherein said 
stiffness parameter unit comprises an iterative processing 



44. The system according to claim 41, wherein said 
stiffness parameter unit comprises an outer iterative process- 
ing unit and an inner iterative processing unit. 

45. The system according to claim 43, wherein said 
iterativ e pi i icessing unit determines said stiffness parameters 
using a lirst order perturbation process. 

46. The system according to claim 43, wherein said 
iterative processing unit determines said stiffness parameters 
using a higher order perturbation process. 
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ABSTRACT 



A method for determining alteration of a region of an earth 
formation surrounding an earth borehole, comprising the 
steps of providing a logging device that is moveable through 
the borehole; transmitting sonic energy into the formation 
and receiving, at a plurality of transmitter-to-receiver 
spacings, sonic energy that has traveled through the 

received sonic energy for the plurality of transmitter-to- 
receiver spacings; determining sonic transit times and dif- 



ttial l 



r the 



espective transmiller-lo- 
r spacings; deriving a test statistic from the 
differential transit times; and determining the presence of 
alteration of a region of the formations from the test statistic. 
An associated apparatus for carrying out the method is also 
described. 
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US 6,526, 

1 

SONIC WELL LOGGING FOR ALTERATION 
DETECTION 

RELATED APPLICATION 

The present application contains subject matter that is 5 
related to subject matter in copending U.S. patent applica- 
tion Ser. No. 09/741574, entitled "Sonic Well Logging for 
< 'baraclcri/ing Earth Formations", Attorney Docket Number 
60.1406, incorporated herein by reference, filed of even date lg 
herewith, and assigned to the same assignee as the present 
application. 

FIELD OF THE INVENTION 

This invention relates to investigation of earth formations, 15 
to a method and apparatus for determining properties of 
earth formations using sonic well logging which can char- 
acterize earth formations exhibiting complex acoustic 
behavior, such as can occur in anisotropic and/or inhomo- 
geneous formations, and to a method and apparatus for 20 
determining alteration of a region of formations surrounding 
an earth borehole. 

BACKGROUND OF THE INVENTION 

It is well known that mechanical disturbances can be used - 
to establish acoustic waves in earth formations surrounding 
a borehole, and the properties of these waves can be mea- 
sured to obtain important information about the formations 
through which the waves have propagated. Parameters of 
compressional, shear and Stoneley waves, such as their 30 
velocity (or its reciprocal, slowness) in the formation and in 
the borehole, can be indicators of formation characteristics 
that help in evaluation of the location and/or producibility of 
hydrocarbon resources. 

An example of a logging device that has been used to 
obtain and analyze sonic logging measurements of forma- 
tions surrounding an earth borehole is called a Dipole Shear 
Sonic Imager ("DSP" — trademark of Schlumberger), and is 
of the general type described in Harrison et al., "Acquisition 
and Analysis of Sonic Waveforms From a Borehole Mono- 
pole And Dipole Source For The Determination Of Com- 
pressional And Shear Speeds And Their Relation To Rock 
Mechanical Properties And Surface Seismic Data", Society 
of Petroleum Engineers, SPE 20557, 1990. In conventional 4J 
use of the DSI logging tool, one can present compressional 
slowness, At c , shear slowness, At,, and Stoneley slowness, 
At st , each as a function of depth, z. [Slowness is the 
reciprocal of velocity and corresponds to the interval transit 
time typically measured by sonic logging tools.] Typically, 
the subsurface formations are considered to be homoge- 
neous and isotropic material, where the compressional and 
shear velocities, V c . and V,, of the formations are only a 
function of depth. 

It is known that formations can be anisotropic, where the 55 
compressional and shear slownesses are a function of 
azimuth, 6. Anisotropy can occur, for example because of 
layered shales, aligned fractures or differences in the mag- 
nitudes of the principal stresses in the formations. 

It is also known that formations may be inhomogeneous, 60 
where the slownesses become a function of radial distance 
r, from the borehole. Inhomogeneity can be caused, for 
example, by mud-shale interactions or by mechanical dam- 
age due to stress concentrations. 

It is among the objectives of the present invention to 65 
provide an improved technique for characterizing earth 
formations exhibiting complex acoustic behavior. It is 
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among the further objects of the invention to provide an 
improved technique and apparatus for detecting alteration of 
a region of earth formations surrounding a borehole. 
SUMMARY OF THE INVENTION 

A form hereof provides a technique for jointly processing 
monopole and dipole measurements from a sonic logging 
device to obtain a more complete characterization of the 
acoustic behavior of formations surrounding an earth bore- 
hole. From an integrated sonic inversion, improved esti- 
mates of slowness, as well as characterization of the near 
wellbore damage region, can be obtained. 

In accordance with a form hereof, a method is set forth for 
determining properties of an earth formation surrounding an 
earth borehole, comprising the following steps: (a) provid- 
ing a logging device that is moveable through the borehole; 
(b) transmitting sonic energy into the formation, receiving, 
at the logging device, sonic energy that has travelled through 
the formation, and producing signals representative of the 
received sonic energy; (c) determining, from the signals, 
whether the formation is anisotropic; (d) determining, from 
the signals, whether said formation is inhomogeneous; and 
(e) responsive to the determinations of steps (c) and (d), 
outputting a characterization of the formation as one of the 
following types: isotropic/homogeneous, anisotropic/ 
homogeneous, isotropic/inhomogeneous, and anisotropic/ 
inhomogeneous. In an embodiment of this form of the 
technique, the step (b) includes transmitting sonic energy 
from a monopole source, and receiving sonic energy from 
the monopole source at a plurality of different transmitter- 
to-receiver spacings on the logging device, and the step (d) 
includes determining whether said formation is inhomoge- 
neous from deviations between signals at different 
transmitter-to-receiver spacings. In this embodiment, the 

step (b) also includes transmitting dipole shear sialic energy, 
receiving dipole shear sonic energy that has travelled 
through the formation, and producing signals representative 
of the received sonic energy over a range of frequencies, and 
the Step (c) includes determining whether the formation is 

anisotropic from said signals. It will be understood that the 
step of transmitting dipole shear sonic energy can comprise 
producing what is commonly referred to as a flexural wave 
by employing a dipole source in the borehole to cause a 
flexing of the borehole wall. 

In a further form hereof, a method is set forth for 
determining properties of an earth formation surrounding an 
earth borehole, comprising the following steps: (a) provid- 
ing a logging device that is moveable through the borehole; 
(b) transmitting monopole and dipole sonic energy from the 

logging device, monopole and dipole sonic energy that has 
travelled through the formations, and producing measure- 
ment signals representative of the received monopole and 
dipole sonic energy; (c) devising a plurality of formation 
models of different complexities; (d) comparing model 
signals from the models with the measurement signals; and 
(e) selecting one of the models based on said comparing 
step. In an embodiment of this form of the invention, the 
plurality of models includes four models, the models, in 
order of increasing complexity, being a homogeneous/ 
isotropic model, a homogeneous/aniso tropic model, an 
inhomogeneous/isotropic model, and an inhomogeneous/ 
anisotropic model. Also in this embodiment, the step (e) 
selection of a mode! lakes into account model complexity as 
well as the results of the comparing step. 

A feature of the new sonic characterization of a form 
hereof is that a joint processing of compressional and shear 



US 6,526,354 B2 



data is being used to indicate the "state" of the formation, in 
the context of inhomogeneity and/or anisotropy. This "state" 
of the formation would be a function of depth. Some of the 
applications fall into two general categories: those requiring 
"undamaged" parameters, and those applications requiring 5 
"damaged" parameters. "Undamaged" parameters (the 
velocities Y p and of the virgin formation) can be used by 
geophysics and petrophysics experts to evaluate the reser- 
voir and overburden rock in the (radii ional manner. With the 
method hereof, there will be improved confidence that the 10 
measurements have a deep enough depth of inv estigation to 
be representative of the undamaged formation. 

Near wellbore ••damaged" parameters are new informa- 
tion of use to drilling and completions engineers. As an 
example, for inhomogeneity in a reservoir caused by stress 15 
concentrations, the location of damaged zones can influence 
the completion and perforating strategy. One can selectively 
perforate to avoid damaged zones. Also, detection of dam- 
age from sonics can have an impact on other logging 
measurements that have shallow depths of investigations, so 20 
suitable adjustments can be made. 

In accordance with a form of the present invention, there 
is set forth a method foi determining alteration of a region 
of an earth formation surrounding an earth borehole, 
prising the following steps: providi 
is moveable through the borehole; 
into the formation and receiving, at a plurality of 
to-receiver spacings on the logging device, sonic energy that 
has traveled through the formation, and producing signals 
representative of the received sonic energy for the plurality 
of transmittcr-to-rccciver spacings: determining, from the 
signals, sonic transit times and differential transit times for 
the respective transmitter-to-receiver spacings; deriving a 
test statistic from the differential transit times; and deter- 
mining the presence of alteration of a region of the forma- 
tions from the test statistic. In a preferred embodiment of 
this form of the invention, the test statistic includes a 
component that depends on the degree to which the differ- 
ential transit times decrease monolonically as a function of 
transmitter-to-reccivcr spacing, and the step of determining 
the presence of alteration of a region of the formations from 
the test statistic comprises comparing the test statistic to a 
threshold. 

An embodiment of the invention employs a test statistic 
T, of the form 



\ "" " 



where D'lT,- arc the individual differential transit limes, DTf 
is the average of the differential transit times. TR, are the 
individual Iransmitler-to-receiver spacings. and m and c arc 



Further features and advantages of the invention will 
become more readily apparent from the following detailed 
description when taken in conjunction with the accompa- 
nying drawings. 

BRIEF DESCRIPTION OF THE DRAWINGS 
FIG. 1 is a diagram, partially in block form, of a type of 



; that t 



i be i 



apparalt 
hereof. 

FIG. 2 is a simplified diagram of a type of downhole 
logging device that can be used in practicing embodiments 



FIG. 3 is a diagram illustrating, i 
placement of hydrophones that can be used at a receiver 
station in the logging device of FIG. 3. 

FIG. 4 is a block diagram of a portion of the electronics 
of the FIG. 2 logging device. 

FIG. 5 is a conceptual diagram showing, on the left, a 
traditional model used in sonic characterization of forma- 
tions and, on the right, a more general model used in a form 
of hereof. 

FIG. 6, which includes FIGS. 6 A and 6B, show simplified 
sonic logging device representations, with a source and 
receivers, and illustrate ray paths from a monopole source in 
a homogeneous formation (FIG. 6A) and in an inhomoge- 
neous formation (FIG. 6B). 

FIG. 7 illustrates a portion of a log of compressional 
slowness at relatively long and relatively short transmitter- 
to-receiver spacings. 

FIG. 8 shows a dipole dispersion curve of slowness as a 
function of frequency. 

FIG. 9, which includes FIGS. 9A, 9B, and 9C, shows 
dipole dispersion curves for situations of formation isotropy 
(FIG. 9A), formation intrinsic anisotropy (FIG. 9B), and 
formation stress-induced anisotropy (FIG. 9C). 

11(1. 10 shows dipole dispersion curves for a model 
,ging device that 25 formation and for measured data, respectively, and illus- 
trates a situation of near wellbore damage. 

FIG. 11, which includes FIG. 11A and FIG. 11B, shows, 
in FIG. 11A, monopole compressional slowness as a func- 
tion of depth for relatively short and relatively long 
transmitter-to-receiver spacings for a case of a homogeneous 
and isotropic formation, and, in FIG. 11B, dipole dispersion 
curves and a model dipole dispersion curve of slowness as 
a function of frequency for the case of a homogeneous and 
isotropic formation. 

FIG. 12, which includes FIG. 12A and FIG. 12B, shows, 
in FIG. 12A, monopole compressional slowness as a func- 
tion of depth for relatively short and relatively long 
transmitter-to-receiver spacings for a case of an inhomogc- 
neous and isotropic formation, and, in FIG. 12B, dipole 
dispersion curves and a model dipole dispersion curve of 
slowness as a function of frequency for the case of an 
inhomogeneous and isotropic formation. 

FIG. 13, which includes FIG. 13A and FIG. 13B, shows, 
in FIG. 13A, monopole compressional slowness as a func- 
tion of depth for relatively short and relatively long 
transmitter-to-rcceiver spacings for a case of a homogeneous 
and anisotropic formation, and, in FIG. 13B, dipole disper- 
sion curves and a model dipole dispersion curve of slowness 
as a function of frequency for the case of a homogeneous 
and anisotropic formation. 

FIG. 14, which includes FIG. 14A and FIG. 14B, shows, 
in FIG. 14A, monopole compressional slowness as a func- 
tion of depth for relatively short and relatively long 
transmitter-to-receiver spacings for a case of an inhomoge- 
neous and anisotropic formation, and, in FIG. 14B, dipole 
dispersion curves and a model dipole dispersion curve of 
slowness as a function of frequency for the case of an 
inhomogeneous and anisotropic formation. 

FIG. 15 shows an example of the model selection, with 
comparison of simple and complex models. 

FIG. 16 illustrates an example of determining the number 
of dispersion curves resulting from dipole data, using a 
likelihood factor approach. 

FIG. 17 is a flow diagram of a routine for programming 
a processor to implement a routine in accordance with an 
embodiment hereof. 
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FIG. 18 is a diagram of an example of a linear radial 
profile of velocity around the borehole. 

FIG. 19 illustrates the variation of transit time with 
transmitter-to-receiver spacing for the linear radial profile of 
FIG. 18. 5 

FIG. 20 shows differential transit time as a function of 
transmitter-to-receiver spacing for the transit times of FIG. 
19. 

FIG. 21 shows best fit outcomes for exemplary data 1Q 
showing how two test statistics can be fitted to data. 

FIG. 22 is a flow diagram of a routine, in accordance with 
an embodiment of the invention, for determining alteration 
around the borehole and virgin formation slowness. 

DETAILED DESCRIPTION 15 
Referring to FIG. 1, there is shown a type of apparatus 
which can be used in practicing embodiments of the inven- 
tion. Subsurface formations 231 are traversed by a borehole 
232 which is typically, although not necessarily, filled with 20 
drilling lluid or mud. A lagging tool 210 is suspended on an 
armored cable 212 and may have optional centralizers (not 
shown). The cable 212 extends up the borehole, over a 
sheave wheel 220 on a derrick 221 to a winch forming part 
of surface equipment 250. Known depth ganging apparatus :5 
(not shown) is provided to measure cable displacement over 
the sheave wheel 220 and accordingly the depth of the 
logging tool 210 111 the borehole 232. A device of a type well 
known in the art is included in the tool 210 to produce a 
signal indicative of orientation of the body of the tool 210. 30 
Processing and interface circuitry within the tool 210 
amplifies, samples and digitizes the tool's information sig- 

eqiiipnicnt 250 via the cable 212. Electrical power and 

generated by the surface equipment 250 and communicated 
via the cable 212 to circuitry provided within the tool 210. 

The surface equipment includes processor subsystem 270 
(which can typically include a microprocessor, memory, 
clock and timing, and input/output functions — not sepa- 40 
rately shown), standard peripheral equipment (not separately 
shown), and recorder 226. 

The logging device 210 may be, for example, of a type 
known as a Dipole Shear Sonic Imager ("DSI" — trademark 
of Schlumberger) generally described in Harrison et al., 45 
"Acquisition and Analysis of Sonic Waveforms From a 
Borehole Monopole and Dipole Source for the Determina- 
tion of Compressional and Shear Speeds and Their Relation 
t ) Rock Mechanical Properties and Surface Seismic Data". 
Society of Petroleum Engineers, SPE 20557, 1990. It will be 5U 
understood, however, that any suitable logging device can he- 
utilized. Further details of the logging device 210 of this 
example are shown in FIG. 2. The logging device 210 
includes crossed dipole transmitters 315 and 320 (only one 
end of dipole 320 being visible) and a monopole transmitter 55 
325, so that waves including compressional, shear, Stoneley, 
and ffexural can be excited. Eight, or other suitable number, 
of spaced apart receiver stations, designated 331 through 
338 each comprise four receiver hydrophones mounted 
azimuthally at ninety degree intervals in the surface of the 63 
cylindrical logging device. FIG. 3 shows the hydrophones, 
designated A, B, C, and D. In an example shown in FIG. 4, 
an X component can be obtained by subtracting the signals 
received at A and C (i.e., A-C), and a Y component can be 
obtained by subtracting the signals received at B and D (i.e., 65 
B-D). With four receiver elements at each receiver station, 
there are a total of thirty two receiver elements in this 
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example. The receiver stations are also configurable for 
monopole reception. 

The transmitter electronics contain a power amplifier and 
switching circuitry capable of driving the two crossed-dipole 
transmitter elements and the monopole element from a 
programmable waveform. Separate waveforms with appro- 
priate shape and frequency content can be used for dipole, 

electronics processes the signals from the 32 individual 
receiver elements located at the eight receiver stations which 
are spaced six inches apart. At each station, four receivers 
are mounted as shown in FIG. 3 which allows measurement 
of the dipole and crossed-dipole waveforms by differencing 
the outputs from opposite receivers, as previously described. 
Summing the outputs of the receivers can be used to produce 
a monopole equivalent signal. As further described in Har- 
rison et al., supra, the receiver electronics multiplexers, 
filters, amplifies and channels the signals from the 32 
receiver elements to N parallel signal paths. These eight 
parallel analog signals are passed to an acquisition electron- 
ics cartridge where eight 12-bil analog-to-digital converters 
digitize the signals from the receiver electronics. The telem- 
etry circuitry passes the digitized information to the earth's 
surface. 

FIG. 4 shows an example of the acquisition signal path in 
block diagram form for one of the eight (or other suitable 
number of) receiver stations, as described in Harrison et al, 
supra. Each receiver has its own charge preamplifier 
(represented at 505). The output of the receivers, odd 
numbered pairs being in-line with the upper dipole trans- 
mitter and even numbered pairs with the lower dipole 
transmitter, passes into both a summing circuit (for mono- 
pole measurements) and a differencing circuit (for dipole 
measurements), as represented at 510. Under software con- 
trol the sum or difference is selected by a multiplexer stage 
(block 520) and the signal passed to one of eight program- 
mable gain amplifier stages (540) and filters (545). The other 
similar channels are represented by block 550. The eight 
parallel analog signals arc passed to eight parallel 12-bit A/D 
converters (represented at 560) where simultaneous wave- 
form digitization is performed. If desired, more bits can, of 
course, be used to advantage. After digitization, the eight 
waveforms pass to the memory section associated with 
downhole processor 580. The processor also provides con- 
trol signals and waveforms to transmitter and receiver 
electronics. An alternate path directs the eight analog 
receiver signals into threshold crossing detection circuitry or 
digital first motion detection, as represented at block 565. 
This circuitry detects the time of all up or down going 
threshold crossings. The digitized waveform data and the 
threshold crossing time data are passed to the surface using 
telemetry' circuitry 590, It will he understood that more 
advanced tool implementations, having further transmitters, 
receivers, and/or transmittcr-to-rcccivcr (T/R) spacings, and 
more powerful processing capabilities, can be used even 
more advantageously, consistent with the principles hereof, 
in practicing embodiments of the invention. 

In the FIG. 2 embodiment, the processing of signals 
recorded uphole can be implemented using a processor 270, 
such as a suitably programmed general purpose digital 
processor with memory and peripherals conventionally pro- 
vided. It will be understood, however, that the processing 
need not be performed at the wellsite, and that signals 
derived at the wellsite can be processed at a remote location. 
It will also be understood that other suitable logging tools 
can be employed in practicing the invention. 

The traditional sonic log presents compressional 
slowness, At c , and shear slowness, At s , and Stoneley 
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slowness, At Jt , as a function of depth, z. This description has 
usually been based on the idea that rock can be described as 
a homogeneous and isotropic material where the compres- 
sional and shear velocities, Vc and Vs, of the formation are 
only a function of depth. It has become recognized however, : 
that this is not a complete characterization of the formation. 
Formations may be anisotropic where the compressional and 
shear slownesses are a function of azimuth, 6. Anisotropy 
can occur because of layered shales, aligned fractures or 
differences in the magnitudes of the principal stresses, t 
Formations may also be inhomogeneous where the slow- 
nesses become a function of radial distance from the 
borehole, r. Inhomogeneity can be caused, for example, by 
mud-shale interactions or by mechanical damage due to 



A form hereof addresses the compressional and shear 
velocity around a borehole, Vc,s(r,0), at a single depth and 
a processing chain to exploit the model. Conceptually, a 
model of velocity around a borehole model is shown in FIG. 
5. The velocities Vc,s(r,9) associated with the model are 
described as only a function of radius and azimuth. A typical 
old model, shown on the left, depicts a homogeneous and 
isotropic formation with Vc,s(r,6)=constant. In the old 
model, the rock around the borehole is assumed to be the 
same as if no borehole 

model hereof, on the right, assumes that the drilling process 
and related stress concentrations and/or mud system may 
have damaged or altered the rock from its pre -drilled state. 
At this stage, consider the general case where the altered, 
changed or damaged rock (see Plona, T. J., Sinha, B. K., 
Winkler, K. W, and D'Angelo, R., "Measurement of Stress 
Direction and Mechanical Damage Around Stressed Bore- 
holes Using Dipole and Microsonic Techniques", SPE 
47234 presented at Eurock '98 (1998)) can have an elliptical 



small graph inset indicates that the sound speed is substan- 
tially constant with radius. This is the traditional concept of 
monopole sonic logging where the assumption is generally 
that the formation is homogeneous. However, when there is 
a gradient of sound speed vs ratlins (i.e., inhomogeneity), it 
is known that sound propagates in curved ray paths. For the 
case indicated in FIG. 6B, where the sound speed increases 
vs radius from the borehole, the ray paths are curved such 
that energy returns to the borehole. In this case, the ray path 
j that arrives at the nearest receiver senses a velocity that is 
indicative of the formation near the wellbore (that is, a 
shallow depth of investigation). In contrast, the ray path that 
arrives at the farthest receiver senses a velocity which in 
indicative of the formation relatively far away from the 
5 borehole (i.e., a relatively deep depth of investigation, not 
substantially affected or perturbed by the borehole). Thus, 
the measured velocity depends on the T/R spacing. The 
difference in the slowness measurement at short and long 
T/R spacings provides an indication of formation inhomo- 

Consider the example of logs of slowness data as in FIG. 
7. The track 710 represents the slowness at the long T/R 
spacing and is generally indicative of the unperturbed or 
undamaged rock. The track 720 represents the slowness at 
dnrdibn^hrbcn"cLTto 25 the short T/R spacing and is characteristically a higher 
imes that the drilling nrocess slowness (i.e. lower velocity). The difference in the two 



curves, at depth levels between about 100 and 150 ft. i 
example, indicates radial variations of velocity (i.e., 
inhomogeneity). 
3 For the dipole shear measurements, the data is acquired in 
the crossed dipole mode, and the data can be used to 
determine the isotropy or anisotropy of the formation. The 
four-component rotation method of Alford can be used. (See 
Alford, R. M., 1986, "Shear Data In the Presence Of 
region around the borehole that has' been perturbed The 35 Azimuthal Anisotropy": 56th Ann. Internal., Soc. Expl. 
elliptical region implies that the formation can be both Geophys., Expanded Abstracts, 476-479)). Slowness is 
anisotropic and inhomogeneous, i.e., Vc,s-V(r,e). evaluated as a function of frequency (i.e., the dispersion 

_ ri1 , , , , . , , , . , , curvet In FIG X, a basic dipole llcxural mode dispersion 

With the rock around the borehole modeled to have both . , . „.„ r f , ,. „ , * f tU 

• r , , . _■• curve is shown, where d refers to the diameter ol the 

radial and a/imuthal variations ol sound speeds, an embodi- 
ment hereof uses both monopole and dipole data from a 4 
sonic tool or tools, processed jointly, to yield a more 
complete characterization of the formation acoustic proper- 
ties. Historically, sonic logging has evaluated compressional 
and shear data independently. In a form of the present 
il and shear data 



joint manner in order to get a more complete ch 
lion of the formation Inputs lo this charactci i/ation 
monopole compressional slowness and the dipok 
slowness. 

For the monopole input, in this form of the 
monopole compressional data is obtained at both relatively 
short transmitter receiver (T R) spacings (for example, 3 to 
5 ft.) and relatively long T R spacings (for example. 8 to 12 
ft.). The data can be obtained using any suitable technique, 
for example with the type of equipment described in con- 
junction with FIGS. 1-4. Reference can be made, for 
example, to the data obtained in Hornby, B. E. "Tomo- 
graphic Reconstruction of Near-Borehole Slowness Using 
Refracted Borehole Sonic Arrivals", Geophysics, 58, 
1726-1738, (1993). 

In the simplified diagrams of FIGS. 6Aand 6B, formation 
is represented at 601, a borehole is represented at 605, a 
logging device is represented at 620, a monopole source is 
represented at 625, and receivers are represented at 628. 

In FIG. 6Athe sound rays are shown as propagating from 
the monopole source to the receivers in straight lines. The 



borehole. There are several points to observe from this 
Figure. First, the dispersion curve approaches the formation 
shear slowness at low frequencies. Second, at low frequen- 
cies (and long wavelength) the dipole signal is sensing 2 to 
3 borehole diameters into the formation. Third, at high 
45 frequencies (and short wavelengths), the dipole signal is 
sensing approximately Vi the borehole diameter into the 
formation. Thus, the information at different frequencies is 
sensing sound speeds at different depths of investigation. 
(See Sinha, B. K., Norris, A. N., and Chang, S. K., 1994, 
"Borehole Flexural Modes In Anisotropic Formations": 
Geophysics, 59, 1037-1052.) 

In FIG. 9, which includes graphs 9A, 9B and 9C, it is 
shown how the dispersion curves indicate the isotropy or 
anisotropy of the formation. In the plot 9A, there is only one 
dispersion curve visible. This is because the two curves from 
the respectively orthogonal dipole pairs (see e.g. FIG. 3) 
substantially overlap. This is the case when the formation is 
isotropic. The dispersion curve approaches the shear speed 
(dashed line) at low frequencies. In the plots 9B and 9C there 
are two dispersion curves, and this indicates that the forma- 
tion is anisotropic. The plot 9B shows the case when the 
formation anisotropy is "intrinsic" and the two dispersion 
curves do not cross. (See Sinha, B. K., Norris, A. N., and 
Chang, S. K., 1994, supra.) The plot 9C shows the case when 
there is stress-induced anisotropy and the two dispersion 
curves cross. (See Sinha, B. K., and Kostek, S., 1996, 
"Stress Induced Azimuthal Anisotropy In Borehole Flexural 
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Waves": Geophysics, 61, no. 6, 1899-1907.) The crossing 
nature ol the stress-induced anisotropy dispersion curves is 
a feature that permits a discrimination to be made between 
the two types of anisi<lr' >py. These stress-induced anisotropy 
effects have been observed in field data. (Sinha B. K., Kane, : 
M., and Frignet, B., "Dipole Dispersion Crossover And 
Sonic Logs In Limestone Reservoirs", Geophysics, 65, No. 
2 (March-April 2000), pp. 390-407). 

In addition to the distinction between isotropy and 



these examples the monopole and dipole data have been 
derived from different wells, the data has been combined in 
this set of examples to illustrate the principle hereof. It will 
be understood that, in practice, the measurements can be 
taken by a single tool (or, less preferably, plural tools in a 
single logging run, or plural logging runs). 

FIGS. 11A and 11B illustrate the homogeneous and iso- 
tropic case. From the monopole data of FIG. 11A, it is seen 
that the compressional slowness from the short T/R n 



anisotropy, disper: 



homogeneity/inhomogeneity. This is illustrated in conjunc- 
tion with FIG. 10. In this Figure, the model data is repre- 
sented at curve 1021 and the measured data is represented at 
1031. The model data can be produced, for example, from 
measured compressional and shear velocities, formation 
mass density, mud density, mud compressional velocity, and 
borehole diameter (see Sinha, B. K., Norris, A. N., Norris, 
A. N, and Chang, S. K., 1994, Borehole Flexural Modes In 
Anisotropic Formations, Geophysics, 59, 1037-1052). If the 
measured data superimposes with the model data, it can be 
concluded that the formation is homogeneous. When the 
measured data deviates at high frequency (as in the Figure), 
it can be concluded that the formation is inhomogeneous. 
Since the deviation occurs at high frequency that corre- 
sponds to probing near to the borehole, this deviation 25 
indicates that there is lower sound speed near the borehole 
wall, i.e., inhomogeneity or damage. 

With the basic realization that the rock around the bore- 
hole may have both radial and azimuthal variations of sound 
speeds, both monopole and dipole data from sonic tools can 30 
be utilized to give a more complete characterization of the 
formations acoustic properties. In a preferred embodiment 
hereof, the sonic characterization includes four categories as 
follows: 

Homogeneous and Isotropic 
Inhomogeneous and Isotropic 
Homogeneous and Anisotropic 
s and Anisotropic 



also yield an indication of 13 surement (the log 1110) substantially corresponds 



i the 

slowness from the long T/R measurement (the log 1120). 
This is one indication of homogeneity. From the dipole data 
of FIG. 11B, two things can be observed. First, the curves 
1151, 1152 (measured data) are from the crossed dipole 
15 acquisition from orthogonal transmitters receivers. Since the 
two dispersion curves lie substantially on top of each other, 
an isotropic formation is indicated. Secondly, the measured 
data is similar to the modeled dispersion curve (1153). The 
model curve is for the homogeneous case. The measured 
10 data fitting the model is a second indicator that the formation 
is homogeneous. Thus, with a joint processing of monopole 
and dipole data, it is determined that the formation at the 
depth of these measurements is both homogeneous and 
isotropic. 

FIGS. 12A and 12B illustrate the inhomogeneous and 
isotropic ease, f rom the monopole data, it can be seen that 
the compressional slowness from the short T/R measure- 
ment (the log 1210) does not correspond to the slowness 
from the long T/R measurement (the log 1220). This is one 
indication of inhomogeneity, as described above. From the 
crossed dipole data of FIG. 12B, it can be observed that the 
dispersion curves 1251, 1252 (measured data) lie substan- 
tially on top of each other. This indicates that the formation 
is isotropic. Next, a determination is made as to whether the 



it the 



odeled dl 



(1253). Since, as seen, the measured data does not fit the 
model, there is a second indicator that the formation is 
inhomogeneous. Thus, with a joint processing of monopole 
and dipole data, it is determined that the formation at the 
In the present embodiment, monopole compressional slow- 40 depth of these measurements is inhomogeneous and isotro- 
nesses and dipole shear slownesses are used in determining pic. 

the appropriate characterizations. From the monopole com- FIGS. 13A and 13B illustrate the homogeneous and 

pressional data, a determination is made as to whether or not anisotropic case. From the monopole data, it can be seen that 
the compressional slowness is a function of transmitter-to- the compressional slowness from the short T/R measure- 
receiver (T/R) spacing. If the compressional slowness is not 45 ment (log 1310) substantially corresponds to the slowness 
a function of T/R spacing, then the formation is deemed from the long T/R measurement (log 1320). This indicates 
homogeneous. If the compressional slowness is a function of 
T/R spacing, particularly, if the slowness decreases with 
increasing T/R spacing, the formation is deemed to be 
inhomogeneous. [Reference can also be made to the tech- 5 
nique described in conjunction with FIGS. 18-22.] From the 
crossed dipole shear data, it can be determined if the 
formation is isotropic (one shear slowness) or anisotropic 
(two shear slownesses). In addition, inhomogeneity can be 
determined from the deviation of the measured dispersion J 
curve from the modeled curve. With these determinations, 
an acoustic description of the formal ion can be given as one 
of the above-listed characterizations: homogeneous 
isotropic; inhomogeneous/isotropic; homogeneous/ 
anisotropic; or inhomogeneous/anisotropic. ( 

The following set of Figures, FIGS. 11-14, show 
examples of types of sonic field data that highlight 
homogeneity/inhomogeneity and/or isotropy/anisotropy in 
formations. The monopole data of these examples is 
extracted from the idealized data of FIG. 3. (See Hornby, t 
supra.) The dipole Jala is field data from the sonic logging 
tool described in conjunction with FIGS. 2-4. Although in 



homogeneity, from llie ciosscd dipole dala, 
observed that the dispersion curves 1351, 1352 (measured 
data) do not lie on top of each other. There is a characteristic 
difference in slowness at low frequency. This indicates that 
the formation is anisotropic. Next, the two measured dis- 
persion curves are nearly parallel, or non-crossing. This 
indicates that this is an intrinsic anisotropic formation. Next, 
a determination is made as to whether the measured disper- 
sion curves substantially fit the modeled dispersion curve 
(1353). Since, as seen, the measured data substantially fits 
the model, there is a second indicator that the formation is 
homogeneous. Thus, it is determined that the formation at 
the depth of these measurements is homogeneous and aniso- 
tropic. 

FIGS. 14A and 14B illustrate the inhomogeneous and 
anisotropic case. This is the most general case. From the 
monopole data, it can be seen that the compressional slow- 
ness from the short T/R measurement (log 1410) does not 
correspond to the slowness from the long T/R measurement 
(log 1420). This is one indication of inhomogeneity. From 
the crossed dipole data, it can be observed that the dispersion 
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curves 1451, 1452 do not lie on top of each other. There is 
a characteristic difference in slowness at low frequency that 
indicates that the formation is anisotropic. In addition, the 
measured dispersion curves do not fit the modeled disper- 
sion curve (1453), another indication of inhomogeneity. 5 
Thus, it is determined that the formation at the depth of these 
measurements is inhomogeneous and anisotropic. 

In a form of the present invention, an integrated inversion 
of monopole, dipole and Stoneley data is performed to select 
a suitable rock model. After classifying the formation using 
one of the available rock models, the processing chain 
subsequently computes answer products appropriate to the 
selected model. In an embodiment hereof, the first stage in 
the processing chain is to test all of the recorded sonic data 
for consistency with four different models of the physics as 
first set forth above in the Table: 

Homogeneous/I so tropic 

Homogeneous/Anisotropic 

Inhomogeneous/Tsotropic 

Inhomogeneous/Anisotropic 20 
The assumed models vary in complexity. The homogeneous 
isotropic model is the simplest model whose physics is well 
understood i.e., there is a complete theoretical understanding 
of sound pn ipagation in a homogeneous isotropic h u illation. 
The inhomogeneous/isotropic and the homogeneous/ 25 
anisotropic models are more complex, requiring more physi- 
cal parameters for their description. For each of these 
models, there is not as complete an understanding of sound 
propagation as there is for the homogeneous isotropic case. 
The inhomogeneous/isotropic model is the most complex 30 
model, requiring the most physical parameters for its 
description, and the least well understood of the four mod- 
els. 

The different proposed rock models have varying 
complexity, and have the property that the same data might 35 
be explained equally well by different models. Sound propa- 
gation in the homogeneous/isotropic model can be viewed as 
a special case of sound propagation in the homogeneous 

anisotropic model (or, indeed, the other two models) with 



US 6,526,354 B2 

12 

criterion used for model selection is based on computing 
Bayesian posterior probabilities for the four different model 
types. FIG. 15 illustrates the qualitative basis of this proce- 
dure. A simple model, with a few parameters, describes a 
smaller family of possible datasets than a more complex 
model with more parameters. If the horizontal axis describes 
the universe of possible datasets, then the simpler Model 1 
in TIG. 15 describes a smaller subset of datasets compared 
with the more complex Model 2. Notice that there can be an 
overlap; that is, certain sets of data can be described by both 
Models 1 and 2. If one assumes a priori uniform probability 
density distributions are imposed for datasets described by 
both models, then it follows that on the region of overlap, the 
simpler Model is. a priori, the more probable one 

This forms the basis of the model selection procedure of 
the present embodiment. The data is tested for goodness of 
fit between a simpler model and a more complex model, 
which subsumes the simple model. The more complex 
model is selected not simply if it fits the data better than the 
simple model; it is only selected if the goodness of fit is 
sufficiently better to warrant using the added complexity of 
the more complex model. In this embodiment, two specific 
statistical tests are applied to select one of the four models. 
First, there is a test for homogeneity versus inhomogeneity, 
such as the above-described type of testing for radial varia- 
tions of velocity, which can be applied to both monopole and 
dipole data. Second there is a test for isotropy versus 
anisolropy, such as described above, applied only to dipole 



As above, the test for inhomogeneity on monopole data 
comprises a test for differences in monopole compressional 
slowness near and far from the borehole. Measured transit 
times, or differential transit times, are tested for deviations 
from a linear trend (predicted by a homogeneous isotropic 
model). Also as above, the test for inhomogeneity on dipole 
data comprises a test for deviation of a measured flexural 
dispersion curve from a library of known llexural dispersion 
curves belonging to the homogeneous isotropic model (see 
e.g. FIG. 10). The measured flexural dispersion curve can be 
the output of a known dispersion analysis algorithm (see, for 



in physical parameters set to special values (e.g. the two 40 example, Ekstrom, M. P., 1995, Dispersion Estimation From 



shear speeds set to the same value). Thus, simply selecting 
a model based on finding the model which best fits the data 
may not always be a good criterion on which to base the 
classification. Qualitatively, a more complex model, with 
more degrees of freedom, can always be adjusted to fit the 4 
data better than a simpler model. 

Differing model complexity also highlights a related 
problem. The more complex models require more physical 
parameters for their description, some of which cannot be 
easily inferred from the data. In some cases, more complex 5 
models require extremely complicated inversions, or may be- 
severely underdetermined by the available data. An example 
is the inhomogeneous isotropic rock model, which is 
o address the case of radial variation of slowness 



Borehole Acoustic Arrays Using A Modified Matrix Pencil 
Algorithm: 29th Asilomar Conf. Signals Sys, and Compt, 
Pacific Grove, Calif., October 31). The measured flexural 
dispersion with error bars is tested for whether it belongs to 
the known family of curves corresponding to the homoge- 
neous isotropic model, c 



Foil 



/ing 



test for 



mho 



homogeneity, a test is performed for isotropy versus anisot- 
ropy. In the present embodiment, and as described in the 
examples above, this test is only performed on dipole data. 
The test can be a statistically based test for the presence of 
one flexural shear (isotropic) versus two flexural shear 
waves (anisotropic). This can be implemented using com- 
parison of a simple model (one wave present) versus a more 
away from the borehole. The inhomogeneous isotropic 55 complex model (two waves present), and selection of the 



model is much rj 
anisotropic model, a 

In a form hereof, the technique for model selection 
overcomes the problem of differing model complexity, and e 
the problem of indeterminacy for more complex models. 
Initially, testing is performed for deviations from the 
homogeneous, isotropic model. This phase relies only on the 
complete understanding of the homogeneous isotropic 
model and does not require a complete theoretical descrip- t 
tion of the other three models, save that they all subsume the 
homogeneous isotropic model. In an embodiment hereof, the 



model using a likelihood factor, as illustrated in FIG. 16, 
which shows an example in which model order 2 indicates 
that the estimate of two dispersion curves fits the data best; 
i.e., has the largest "evidence". An alternative is to use 
cross-line energy as a test statistic and employ a threshold to 
decide between anisotropy and isotropy; a technique that is 
already used as a qualitative indicator of anisotropy. In 
addition, if Stoneley wave slowness data is available, that 
information can be incorporated as a constraint on the shear 
slowness. 

FIG. 17 is a tlovv diagram that represents the processing 
chain, and which can be used in programming a suitable 
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processor, such as the processor 270 of the surface equip- 
ment of FIG. 1 or a remote processor, in practicing an 
embodiment hereof. The block 1710 represents the collec- 
tion of the data signals. A downhole processor could also 
perform at least part of the technique. The data may, for 
example, be collected and stored using the logging apparatus 
described in conjunction with FIGS. 1-4, although it will be 
understood that other suitable equipment can be utilized. 
The data in the present embodiment includes both monopole 
and dipole broadband measurements. 

The block 1715 represents the generating of the logs (e.g. 
in FIGS. 7, 11A, 12A, 13A and 14A) and the actual and 
model dispersion curves (e.g. in FIGS. 8-10, 11B, 12B, 13B 
and 14B) from the collected measurement data. The block 
1720 represents determination of the monopole slowness 
difference for the relevant depth range, for example, the 
slowness difference of FIGS. 11 A, 12 A, 13 A or 14A. The 
block 1730 represents performance of classification of the 
dispersion curves as was explained in conjunction with 
FIGS. 11B, 12B, 13B and 14B, and which can utilize, for 
example, the type of technique described in conjunction 
with FIG. 16. The decision block 1740 represents the 
determination of whether the formation is homogeneous. In 
the present embodiment, the prior determination of block 
1720 is utilized, together with a predetermined threshold, to 



In an embodiment of the invention, the problem of 
detecting radial alteration (a radial decrease of slowness) in 
the near wellbore region is recast as a problem of detecting 
a variation in measured sli iwness al different TR spacings. It 
5 can be noted that the approach of this embodiment applies 
only to the case where the slowness decreases with increas- 
ing radius; that is, away from the borehole. 

Consider first an idealized case, with a perfect borehole 
and with no lateral variations in slowness in the formations, 
in FIG. 18 shows a linear radial profile of velocity around the 
borehole. In this example the velocity at the borehole wall 
(r=0.35 ft) is about 9000 ft/sec and the velocity in the virgin 
formation (starting at r=l .5 ft) is about 10.000 ft sec. with 
the velocity increasing linearly in the alteration region, 
is FIG. 19 illustrates (solid line curve 1910) the variation of 
transit time (TT) with TR spacing for the linear radial profile 
of FIG. 18. (The corresponding curves, actually straight 
lines, shown in dashed line, for homogeneous profiles with 
borehole wall velocity (1920) and virgin formation (1930) 
20 are presented for comparison.) As seen, the curve 1910 
deviates from linear moveout with TR spacing. This is 
readily observed in FIG. 20, which shows the differential 
transit time (DTT) as a function of TR spacing for Vz ft. 



spacings. The derivative of the transit time with respect to 

determination of homogeneity inhomogeneily, for 25 TR spacing or measured slowness shows a monotonic 

example by determining whether the difference in slowness decline until about TR C which, in both FIGS. 19 and 20, is 

(as between the relatively short and relatively long T/R the minimum TR spacing at which modeled sonic energv 

spacing measurements of slowness), and deeming the for- rays penetrate the virgin formation. 

mation homogeneous if the difference does not exceed the In actualit y ; the TR spacing is not continuously varied, but 

predetermined threshold. Il will be understood, however. , r; nher is available f (1 r cerlain disciele values w hich depend 

that the degree of inhomogeneity, which will be related to on thc transm it tcr an d receiver placements on the logging 

the magnitude of the stated difference, can also be output. tooL Thcrcl ' orc . as in 1IG . 2 (). differential transit limes 

[Also, as noted above, a dipole dispersion versus model can derived from di ff erences between transit times for adjacent 

also be used as an inhomogeneity indicator.] Determination rccc i VC rs for a given transmittcr, can be plotted against thc 

is then made (decision block 1750. shown as two separate 35 average TR spacing of the receivers (' h 11., in this example), 

blocks, 1750A and 1750B in the diagram of FIG. 17) is then It can be noted that the differential transit time exhibits a 

entered for a determination of whether the formation is clear signature of a monotonically decreasing trend leveling 

isotropic or anisotropic. In the present embodiment, this out t0 a constant value wnich is the virgin s i 0W ness. This 

determination can depend on the previously performed provides a basis for devising a detector for radial alteration, 

classification of dispersion curves (block 1730) using, tor 40 , n ^ , embodiment an alteration deteclor is estab . 

example, the rules described in conjunction with the ^ tQ soK . c ;| hvpolhtMS lcslini , proMcm bascd on lh e 

examples of FIGS^ 11B, 12B, 13B, and 14B. It will be above . described differential transit times. It seeks to decide 

understood that the degree of amsotropy can also be be(ween ^ folk)wi two hypotheses , H0 H1: 
i] u a ill died, based, lor example, on the extent ol sepaialion ol 

the dispersion curves. As seen in FIG. 17, depending on the 45 H0: DTT— constant 

outputs of decision blocks 1750A and 1750B, the charac- HI: DTT — decreasing trend with TR. 

leri/.alion of the formation is made as one of: homogeneous In this embodiment there is no assumption of a particular 

isotropic, homogeneous/anisotropic, inhomogeneous/ radial slowness profile, only that slowness decreases with 

isotropic, or inhomogeneous/anisotropic. These increasing radius so as to produce a decreasing trend in 

characterizations can be read out, and the next depth level of 50 measured slowness with TR spacing. The detected signals 

interest processed in similar manner. include measurement errors, and a noise component is 

In accordance with an embodiment of the invention, included in the model, 

detection of inhomogeneity (e.g. detection of alteration To a first order, a decreasing trend can be fitted by a line 



around the borehole resulting from fluid invasion or from 
stress concentration around the wellbore) includes determi- 5 
nation of changes in the measured slowness as the 
lransmiller-lo-reevivcr (TR) spacing is varied. An existing 
model would naturally predict, for a ramp-shaped velocity 
profile, a decrease in measured slowness as the TR spacing 
is increased. However, Applicant has noted that a problem 6 
with using transit time estimates to detect changes in slow- 
ness is complicated by the presence of borehole and forma- 
tion effects as well as gross errors that can occur in transit 
time estimates as a result, for example, of cycle skips. A 
feature hereof is the robust determination of alteration 6 
around the borehole, even in the presence of these compli- 
cating effects. 



with a negative slope. For t 
can be cast as follows: 
HO: DTT=s+N 
HI: DTT=mTR+c+N, m<0 
Where s, c are constants, m is 
measurement error (noise). 

Although the foregoing does r 
the actual DTT relationship i 



the behav 



;, the above hypothes 



the slope, and N is the 



ily precisely fit 
alteration, it does capture 



j first order and is relatively easy to solve for 



The Generalized Likelihood Ratio Test (GLRT) for the 
above problem with independent identically distributed 
(i.i.d.) normal errors yields the following test statistic, T v 



15 
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(5) 



If the errors are not i.i.d, the above statistic can be modified 
by suitably incorporating the covariance matrix of the errors 
as weighting coefficients in the terms of the sums, i.e., 

Yj.WiXDTT, -dtT) 2 (1A) 
'' nii]i);,M' ii (D7T i -»i™ i -c) 2 

Where W=R _1 is the inverse of the covariance matrix R, W,-, 
being the weighting coefficients on the diagonal of W. 

The test consists of comparing the statistic to a threshold 
t which is chosen based on a desired level of false alarms 
(detection of altered zone where there is none) and choosing 
HO when the threshold is exceeded. This can be written as 

H, (2) 
T, r which is interpreted as T, < t => H 0 , Ti > r => 



In a further form of the present embodiment of the invention, 
the test statistic can be refined further by noting that the next 
order fit, especially when there are constant regions in the 
DTT curve, is a piecewise fit with a monotonic decreasing 
trend in the smaller TR spacings and with a constant for the 
longer TR spacings, such that the whole function is mono- 
tonic decreasing. This yields the following test statistic 

^ ; (077;-7J7T) 2 ( 3 ) 

T ' = r c n JSio^^n Si {DTIi ~ (77?i))2 



^ w _|mx + Cl ,x<* c (4) 

III practice, R, is chosen to fall in the middle third of the 
range of TR spacings so as to have enough points for 
statistical stability. FIG. 21 illustrates the way the straight 
line and the piecewise straight line are fitted to the data (the 
star points representing simulated data, including noise) 
while obtaining T 1 and T 2 . The Figure shows best fit 
outcomes for T 1 (curve 2110) and T 2 (curve 2120). The 
expected error free measurement is shown at curve 2130. 

Since there are two test statistics T 1 and T 2 , each of which 
is appropriate in a particular scenario, they can be combined, 
and the problem can be expressed as that of choosing 
between three hypotheses: 

HO: DTT— constant 

Hll: DTT— decreasing linear trend with TR 

H12: DTT — decreasing piecewise linear and constant 

This problem can be solved using the GLRT approach by 
introducing two threshold numbers x 0 and x 12 and imple- 
menting 



Numerical simulations using the present embodiment have 
shown that the performance of the more complex detector 
based on both J 1 and T 2 for detecting alteration may not be 
l , significantly superior to that based on T alone. Accordingly, 
a variation of the present embodiment uses Tj alone to detect 
alteration and uses the ratio of T 2 to T-, to detect the virgin 

The threshold numbers x 0 and x 12 are preferably chosen to 
control the probabilities of false alarm (detecting alteration 

>r ' where none is present) and misclassilieation (finding a 
leveling out when the actual trend is a continuous decrease). 
Thus, it is advantageous to know the distributions of these 
test statistics. It is possible to characterize the distribution of 
the test statistic T 1 under HO with Gaussian errors. It can be 

20 shown to have a distribution which is a (V2, Vi) mixture of a 
point mass at 1 and a shifted F-distribution. This does not 
depend on the variance of the error and therefore allows a 
pre -computation of the threshold for fixing the false alarm 
rate. The distribution of the test statistic T 2 is more complex 

25 and, at present, can be selected on the basis of tests and 
heuristics. 

It is seen from equation (2) that the best fit for the T 2 
statistic incorporates a constant fit at long TR spacings. This 
provides an estimate of the virgin slowness when the test for 

J0 virgin zone picks H12. Referring to FIG. 22, there is shown 
a flow diagram of a routine for determination of alteration 
detection that can be used, inter alia, as an input to homo- 
geneity determination (e.g. blocks 1720 and 1740) of FIG. 
17. As above, the block 1710 represents collection of the 
logging data. In the present embodiment, compressional 

35 transit time measurements are derived from several different 
TR spacings. As above noted, such as in conjunction with 
FIG. 4, the data is digitized and processed (block 2215). As 
represented by the block 2225 blocks of data are subjected 
to gross error filtering, which removes data points that are 

40 grossly inconsistent with neighboring data. Wavelet scale 
filtering can then be implemented, as represented by the 
block 2235. Depth alignment, from measurements at differ- 
ent spacings, can then be performed, as represented by the 
block 2250, and DTT versus TR is obtained (e.g. in the 

45 graph of FIG. 19). The constants s and c are selected and the 
negative slope ill is initialized, as represented by the block 
2260. Also, the threshold x, x 0 , and x n: are selected (block 
2265). Then, as represented by the block 2275, the test 
statistics T 1 and/or T 2 can then be computed, in accordance 

50 with equations (1), (1A), and (3), and the hypotheses from 
among HO, HI and/or HO, Hll, and H12 can be chosen in 
accordance with the threshold determinations (2) and or (5). 
This can be used to provide a general and/or specific (e.g. 
from the quantified test statistics, which can also be pro- 

55 duced using modified constants and slopes, as indicated) 
indications of alteration, as an inhomogeneity indicator. 
Also, as seen from equation (2), the best fit for the test 
statistic T 2 incorporates a constant fit at long TR spacings. 
This provides an estimate of the virgin slowness when the 

60 test for the virgin zone selects hypothesis H 12 . Thus, for 
instance, in the example of FIG. 21, the piecewise fit has a 
constant slowness (as a function of TR spacing) in the virgin 
zone of about 101 usec/ft. 

The foregoing technique is illustrated in terms of com- 

65 pressional waves, but can be applicable to other sonic 
waves, an example being shear headwaves in a fast forma- 
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What is claimed is: 

1. A method for determining alteration of a region of an 
earth formation surrounding an earth borehole, comprising 
the steps of: 

providing a logging device that is moveable through the : 
borehole; 



i the 



min (D7T; -mTR, - 



where DTT, are the individual differential transit times, 
DTT is the average of the differential transit times, TR,- < 
are the individual transmitter-to-receiver spacings, and 
m and c are constants. 



YjADTTi-DTTf 



form alio 

receiving, at a plurality of tr 

ings on said logging device, sonic energy that 
traveled through the formation, and producing signals i 
representative of the received sonic energy for said 
plurality of transmitter-to-receiver spacings; 

determining, from said signals, sonic transit limes and 
differential transit times for the respective transmitter- 
to-receiver spacings; 15 

deriving a test statistic from said differential transit times; 

determining the presence of alteration of a region of the 
formations from said test statistic. 

2. The method as defined by claim 1, wherein said 20 
determined sonic transit times and differential transit times 
are sonic compressional transit times and sonic compres- 
sional differential transit times. 

3. The method as defined by claim 1, wherein said step of 
determining the presence of alteration of a region of the 25 
formal ions from said test statistic comprises comparing said 
test statistic to a threshold. 

4. The method as defined by claim 2, wherein said step of 
determining the presence of alteration of a region of the 
formations from said test statistic comprises comparing said 30 

5. The method as defined by claim 1, wherein said test 
statistic includes a component that depends on the degree to 
which the differential transit times decrease monotonically 
as a function of transmitter-to-receiver spacing. 

6. The method as defined by claim 2, wherein said test 
statistic includes a component that depends on the degree to 
which the differential transit times decrease monotonically 
as a function of transmitter-to-receiver spacing. 

7. The method as defined by claim 2, wherein said test 40 
statistic includes a component that depends on the degree to 
which the differential transit times, as a function of 
transmitter-to-receiver spacing, corresponds to a line of 
negative slope. 

8. The method as defined by claim 2, wherein said test 4 
statistic includes components that depend on the degree to 
which the differential transit times, as a function of 
transmitter-to-receiver spacing, corresponds to a line of 
negative slope followed by a line of substantially zero slope. 

9. The method as defined by claim 2, further comprising 5 
determining the compressional slowness of the virgin earth 
formation from said test statistic. 

10. The method as defined by claim 8, further comprising 
determining the compressional slowness of the virgin earth 
formation from said test statistic. 5 

11. The method as defined by claim 1, wherein said test 
statistic T, is of the form 



where DTT, are the individual differential transit times, 
DTT is the average of the differential transit limes, TR, 
are the individual transmitter-to-receiver spacings, and 
m and c are constants. 

13. The method as defined by claim 2, wherein said test 
statistic Tj is of the form 



Y J .^iADTT i -Dnf 

II ..I/),'/. -;,IK. - 



w here D TT, are the individual differential transit times, 
DTT is the average of the differential transit times, TR,. 
are the individual transmitter-to-receivcr spacings, m 
and c are constants, and W„ are weighing coefficients 
on the diagonal of the inverse of the covariance matrix. 

14. The method as defined by claim 2, wherein said test 
statistic Tj is of the form 



(DTT, - DTT) 2 



where DTT, are the individual differential transit times, 
DTT is the average of the differential transit times, TR, 
are the individual transmitter-to-receiver spacings, and 
m, c 0 , Cj, and R c are constants. 

15. Apparatus for determining alteration of a region of an 
earth formation surrounding an earth borehole, comprising: 

a logging device that is moveable through the borehole; 

means on said logging device for transmitting sonic 
energy into the formation and receiving, at a plurality 
of transmitter-to-receiver spacings on said logging 
device, sonic energy that has traveled through the 
formation, and for producing signals representative of 
the received sonic energy for said plurality of 
transmitter-to-receiver spacings; 

means for determining, from said signals, sonic transit 
times and differential transit times for the respective 
transmitter-to-receiver spacings; 

means for deriving a test statistic from said differential 
transit times; and 

means for determining the presence of alteration of a 
region of the formations from said test statistic. 

16. Apparatus as defined by e Iaim 15, wherein said means 
for determining the presence of alteration of a region of the 
formations from said test statistic comprises means for 
comparing said lest statistic to a threshold. 

17. A method for determining alteration of a region of the 
earth formation, for use in conjunction with a technique for 
sonic logging of an earth formation that includes: providing 
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a logging device that is moveable through the borehole; 
transmitting sonic energy into the formation and receiving, 
at a plurality of transmitter-to-receiver spacings, sonic 
energy that has traveled through the formation, and produc- 
ing signals representative of the received sonic energy for : 
said plurality of transmitter-tu-receiver spacings; compris- 
ing the steps of: 

determining, from said signals, sonic transit times and 
differential transit times for the respective transmitter- ^ 
to-receiver spacings, 

deriving a test statistic from said differential transit times; 

determining the pre: 
formations from ! 
18. The method as 
determined si 



(DTT; -DTT) 2 



Y,s DTT >- 



:e of alteration of a region of the 
test statistic. 

ined by claim 17, wherein said 
mes and differential transit times 
ransit times and sonic compres- 



c compressional t 
sional differential transit times. 

19. The method as defined by claim 17, wherein said step 20 
of determining the presence of alteration of a region of the 
formations from said test statistic comprises comparing said 
lest statistic to a threshold. 

20. The method as defined by claim 17, wherein said test 
statistic includes a component that depends on the degree to 25 
which the differential transit times decrease monotonically 

21. The method as defined by claim 17, wherein said test 
statistic includes a component thai depends on the degree lo 
which the differential transit times, as a function of 30 
lransmiller-lo-reccivcr spacing, corresponds to a line of 
negative slope. 

22. The method as defined by claim 17, wherein said test 
statistic includes components that depend on the degree to 
which the differential transit times, as a function of 35 
transmitter-to-receiver spacing, corresponds to a line of 
negative slope followed by a line of substantially zero slope. 

23. The method as defined by claim 18, further compris- 
ing determining the compressional slowness of the virgin 
earth formation from said test statistic. 40 

24. The method as defined by claim 17, wherein said test 
statistic Tj is of the form 



where DTT,- are the individual differential transit i 

DTT is the average of the differential transit times, TR, su ne g at j ve slope. 

are the individual transmitter-to-receiver spacings, and 

m and c are constants. 
25. The method as defined by claim 17, wherein said test 
statistic T 2 is of the form 



Z; (DTT;-P R _ , :/A'. if 



where DTT,. are the individual differential transit times. 
DTT is the average of the differential transit times, TR, 
are the individual transmitter-to-receiver spacings, and 
m, c 0 , c 1; and R c are constants. 

27. A method for determining whether a region of an earth 
formation surrounding an earth borehole is homogeneous, 
comprising the steps of: 

providing a logging device that is moveable through the 
borehole, 

transmitting sonic energy into the formation and 
receiving, at a plurality of transmitter-to-receiver spac- 
ings on saiel logging device, sonic energy thai has 
traveled through the formation, and producing signals 
representative of the received sonic energy for said 
plurality of transmitter-to-receiver spacings; 

determining, from said signals, sonic transit times and 
differential transit times for the respective transmitter- 
to-rcccivcr spacings; 

deriving a test statistic from said differential transit times; 

determining, from said test statistic, whether said region 
of the formation is homogeneous. 

28. The method as defined by claim 27, wherein said 
determined sonic transit times and differential transit times 
are sonic compressional transit times and sonic compres- 
sional differential transit times. 

29. The method as defined by claim 27, wherein said step 
of determining, from said test statistic, whether said region 
of the formation is homogeneous, comprises comparing said 
test statistic to a threshold. 

30. The method as defined by claim 27, wherein said test 
statistic includes a component that depends on the degree to 
which the differential transit times decrease monotonically 
as a function of transmitter-to-receiver spacing. 

31. The method as defined by claim 27, wherein said test 
statistic includes a component that depends on the degree to 
which the differential transit times, as a function of 
transmitter-to-receiver spacing, corresponds to a line of 



Yj.WiXDTTi -DTT) 2 
nin X It i /)/■/■ -»;/'/>' -,' ; 



32. The method as defined by claim 27, wherein said test 
statistic includes components that depend on the degree to 
which the differential transit times, as a function of 
transmitter-to-receiver spacing, corresponds to a line of 
negative slope followed by a line of substantially zero slope. 

33. The method as defined by claim 27, wherein said test 
statistic Tj is of the form 



where DTT,- arc the individual differential transit times, 
DTT is the average of the differential transit times, TR, : 
are the individual transmiller-to-receiver spacings. m 
and c are constants, and \V„ are weighing coefficients 
on the diagonal of the inverse of the covariance matrix, t 

26. The method as defined by claim 17, wherein said test 
statistic Tj is of the form 



Y,. (.DTT, -DTT) 2 
~ mill X. "■)■/•/. mIK, 



where DTT,- are the individual differential transit times, 
DTT is the average of the differential transit times, TR, 
are the individual transmitter-to-receiver spacings, and 
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34. The method as defined by claim 27, wherein said test 
statistic T-l is of the form 
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where DTT ; are the individual differential transit times, 
DTT is the average of the differential transit times, TR f 
are the individual transmitter-to-receiver spacings, m 
and c are constants, and W„ are weighing coefficients 
on the diagonal of the inverse of the covariance matrix. ^ 

35. The method as defined by claim 27, wherein said test 
statistic Tj is of the form 



where DTT, are the individual differential transit times, 
DTT is the average of the differential transit times, TR t - 
are the individual transmitter-to-receiver spacings. and 
m, c 0 , c„ and R r are constants. 



